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ABSTRACT 


CFDSHIP-IOWA  is  a  general-purpose  unsteady  Reynolds-averaged  Navier-Stokes  CFD 
eode  that  has  been  developed,  over  the  past  10  years,  to  handle  a  broad  range  of  ship 
hydrodynamies  problems.  Originally  designed  to  support  both  thesis  and  projeet  researeh  in  the 
areas  of  resistanee  and  propulsion,  it  has  been  sueeessfully  transitioned  to  Navy  and  university 
laboratories  and  industry,  and  has  reeently  been  extended  to  unsteady  applieations  sueh  as 
seakeeping  and  maneuvering.  It  was  developed  following  a  modern  software-development 
philosophy,  whieh  was  based  upon  open  souree,  revision  eontrol,  modular  eoding  using  Fortran 
90/95,  liberal  use  of  eomments,  and  an  easy  to  understand  arehiteeture  whieh  enables  model 
development  by  users. 

Purpose  of  this  report  is  to  provide:  detailed  doeumentation  of  the  modeling,  numerieal 
methods,  and  eode  development;  user  instruetions  on  ereating  input  files  and  post-proeessing; 
reeommended  proeedures  for  verifieation  and  validation;  and  an  example  simulation  for 
CFDSHIP-IOWA  v3.03.  As  a  framework  for  aehieving  sueeessful  simulations,  an  approaeh 
based  upon  formulation  of  an  initial  boundary  value  problem  and  exeeution  of  a  well-defined 
CFD  proeess  is  developed  and  followed  throughout  the  report. 

Example  simulation  and  other  reeent  applieations  demonstrate  the  eapability  of 
CFDSHIP-IOWA  to  simulate  praetieal  ship  hydrodynamies  problems.  Sueeessful  use  in  both 
thesis  and  projeet  researeh  and  transition  to  other  organizations  demonstrates  the  sueeess  of  the 
overall  design  objeetives.  With  inereasing  use  of  CFD  in  design  proeess,  it  is  expeeted  that 
CFDSHIP-IOWA  will  serve  as  a  platform  for  simulation-based  design  and  optimization  of  future 
naval  vehieles. 
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FONT  CONVENTIONS 


The  following  conventions  are  used  in  this  report; 

Italic 

is  used  for  variables  and  mathematical  symbols 

Constant  Width 

is  used  for  programs  and  procedures,  input  variables,  and  in  examples  to  show  the 
contents  of  files  or  the  output  from  commands. 

Constant  Width  Bold 

is  used  in  examples  to  show  commands  or  other  text  that  should  be  typed  literally  by  the 
user.  (For  example,  mpirun  -np  1  cf  d_ship  means  to  type  "mpirun  -np  1 
cfd_ship"  exactly  as  it  appears  in  the  text  or  the  example) 

Constant  Width  Italic 

is  used  in  examples  to  show  variables  for  which  a  context-specific  substitution  should  be 
made.  (The  variable  filename,  for  example,  would  be  replaced  by  some  actual 
filename). 

% 


is  the  UNIX  C  shell  prompt 


1.  INTRODUCTION 

Reynolds-averaged  Navier-Stokes  (RANS)  eomputational  fluid  dynamies  (CFD)  eodes 
have  matured  for  most  diseiplines  and  are  rapidly  being  integrated  into  the  design  process  such 
that  the  reality  of  physics-based  simulation  based  design  seems  imminent.  General-purpose 
research  codes  are  available  for  many  engineering  applications  such  as  aerospace  (Meakin  and 
Wissink,  1999;  Bush  et  al.,  1998),  ship  hydrodynamics  (Hyams  et  al.,  2000;  Larsson  et  al., 
2000),  and  turbomachinery  (Chima,  2001;  Hall  et  ah,  1999)  whereas  other  applications,  such  as 
automobile  and  industrial  processes,  primarily  take  advantage  of  commercial  codes.  Most  codes, 
especially  commercial  ones,  can  handle  more  than  one  application.  Code  development  has 
evolved  from  Ph.D.  thesis  projects  to  dedicated  groups  at  academic  institutes,  government  and 
industry  laboratories,  and  commercial  companies.  These  groups  struggle  to  keep  pace  with  the 
requirements  of  the  design  community  by  addressing  issues  of  modeling,  numerical  methods, 
high  performance  computing  (HPC),  structured  and  unstructured  grids  and  grid  generation,  and 
pre  and  post  processing.  Other  pace  setting  issues  include  lack  of  trained  users  and  consensus  on 
quality  assessment  verification  and  validation  (V&V)  methodology  and  procedures. 

Present  interest  is  in  ship  hydrodynamics,  which  has  unique  features  in  comparison  to 

related  applications  due  large  Reynolds  number  (Re)«  10^;  small  Mach  number  (incompressible 
flow);  tanker,  cargo/container,  and  combatant  geometries;  ballast,  motions,  maneuvering, 
restricted  water,  and  ambient  waves  operating  and  environmental  conditions;  Froude  number  (Fr) 
and  free-surface  effects  (waves,  spray,  breaking,  near-surface  turbulence,  and  boundary-layer 
and  wake  and  vortex  interactions);  and  propulsor-body  interactions  and  cavitation.  Detailed 
physics  vary  considerably  depending  on  geometry  and  operating  and  environmental  conditions. 

The  status  of  ship  hydrodynamics  CFD  for  steady  flow  design  conditions  was  assessed  at 
the  recent  Gothenburg  2000  Workshop  on  CFD  in  Ship  Hydrodynamics  (Larsson  et  ah,  2000). 
Twenty  groups  representing  16  (8  academic  and  8  industrial)  institutes  and  1  commercial  CFD 
code  company  from  1 1  countries  submitted  results  for  one  or  more  of  three  test  cases  for  modern 
tanker,  container,  and  surface-combatant  hull  forms  with  validation  focusing  on,  respectively, 
turbulence  modeling  and  full-scale  Re;  free-surface  effects  and  propeller-hull  interaction;  and 
free-surface  effects.  Verification  was  required  for  each  group  for  at  least  one  of  the  test  cases 
and  the  V&V  approach  of  Stem  et  al.  (2001)  was  recommended.  Most  codes  used  2  or  more 
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equation  turbulence  models  and  had  free-surface  capability;  finite  volume  or  difference  and  2"‘*- 
and  3’^‘*-order  accurate  numerical  methods;  multi  or  single  block  structured  grids;  and  either 
pressure  Poisson-equation  or  artificial  compressibility  formulations.  Few  codes  included 
capability  for  propulsor  modeling,  unstructured  grids,  or  parallel  computing.  Most  groups 
conducted  verification  for  resistance  Ct  and  many  followed  recommended  procedures;  however, 
there  were  some  problems  due  to  solutions  far  from  asymptotic  range  and  lack  of  experience 
with  detailed  verification  procedures,  especially  for  practical  applications.  Seven  codes  showed 
grid  convergence  for  Ct  with  average  number  of  grid  points  1.5M  and  simulation  numerical 
uncertainty  Usn=^-65%.  Some  groups  showed  oscillatory  convergence  and  variability  between 
grid  studies,  which  indicates  need  for  finer  grids.  Quantitative  validation  was  performed  for  Ct. 
The  average  comparison  error  was  E=5%,  which  approximately  equals  the  validation  uncertainty 
[/[^5%  such  that  the  codes  are  approximately  validated  at  the  5%  level,  which  interestingly  is 
also  equal  to  the  coefficient  of  variation  for  Ct.  The  average  experimental  uncertainty  was 
L/d=1.6%.  Only  qualitative  validation  was  performed  for  point  variables.  The  Reynolds  stress 
turbulence  models  performed  best,  however  2  equation  models  were  also  surprisingly  good. 
Both  surface  tracking  and  capturing  methods  showed  good  free  surface  results. 

Advancements  for  off-design  problems  and  unsteady  flow,  although  warranted  since 
previous  CFD  Workshop  Tokyo  (Kodama,  1994),  are  still  relatively  rare.  For  off-design  yaw, 
steady  turn,  and  shallow  water  problems,  steady  flow  methods  can  still  be  used  and  have  shown 
fairly  good  agreement  with  data,  although  issues  remain  as  to  resolution  of  steep  and  breaking 
waves  and  body  and  wave-induced  vortices  (Tahara  et  ah,  1998;  Alessandrini  and  Delhommeau, 
1998;  Flochbaum,  1999;  Di  Mascio  and  Campana,  1999;  Ohmori  et  ah,  2000).  For  unsteady 
flow  problems,  studies  are  very  limited  partly  due  to  the  fact  that  not  all  steady  flow  methods  are 
easily  extended  to  unsteady  flow.  Ohmori  (1998)  performs  simulations  of  unsteady  combined 
sway  and  yaw  motions  as  per  captive  model  testing  in  towing  tanks  using  planar  motion 
mechanisms;  however,  the  free-surface  is  neglected,  i.e.,  simulations  are  for  the  so-called  double 
body  zero  Froude  number  problem.  Beddhu  et  al.  (1998),  Gentaz  et  al.  (1999),  and  Yeung  et  al. 
(2000)  perform  simulations  for  forced  motions,  and  Sato  et  al.  (1999)  and  Rhee  and  Stern  (2001) 
perform  simulations  for  motions  in  regular  head  waves.  With  capability  for  more  practical 
geometries  and  conditions,  coupled  with  continuing  improvements  in  FIPC  resources,  the 
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promise  of  CFD-based  optimization  will  soon  be  realized  (Tahara  et  al.,  2000;  and  Dreyer, 

2002). 

This  report  provides  doeumentation  of  code  development  of  CFDSHIP-IOWA,  which  is  a 
general-purpose  parallel  unsteady  RANS  ship  hydrodynamics  code.  It  is  intended  for  use  both  in 
research  and  design  at  universities,  and  industrial  and  governmental  laboratories.  The  approach 
includes  2-equation  turbulence,  free-surface  tracking,  and  body-force  propulsor  modeling; 
structured  overset-grid,  higher-order  finite-difference,  and  pressure-implicit  split-operator 
(PISO),  numerical  methods;  parallel  and  portable  high  performance  computing  (FIPC);  and  open 
source,  commented,  and  modular  programming  with  revision  control;  and  web  site  distribution 
(http://www.iihr.uiowa.edu/~cfdship).  The  current  version  of  CFDSFIIP-IOWA  has  benefited 
from  predecessor  thesis  codes  (Tahara  and  Stem,  1996;  Rhee  and  Stem  2001),  but  as  a  complete 
package  represents  significantly  improved  modeling,  numerical  methods,  FIPC,  and  overall  code 
development  effort.  Current  version  is  primarily  intended  for  steady  and  unsteady  resistance  and 
propulsion  simulations,  including  option  of  with  or  without  free  surface  and  body  force  modeling 
or  complete  propulsor.  Forthcoming  versions  will  include  capability  for  body  motions  enabling 
steady  and  unsteady  ship  motions  and  maneuvering  simulations. 

Code  development  was  done  over  period  of  last  10  years  during  which  time  earlier 
versions  were  released  (Paterson  et  ah,  1998;  Wilson  et  ah,  1998)  and  used  as  intended  for 
numerous  applications  including:  turbulence  modeling  (Walker,  2000;  Gill,  2000);  two  phase 
flow  modeling  (Larreteguy,  et  ah,  1998);  wave-induced  separation  (Kandysami,  2001); 
maneuvering  (Simonsen  and  Stem,  2003);  surface-ship  motions  (Weymouth  et  ah,  2003;  Wilson 
and  Stern,  2002);  optimization  (Tahara  et  ah,  2000);  propulsor  flows  (Kim  et  ah,  2003;  Chen, 
2000);  and  preliminary  industrial  design  of  DD21,  i.e.,  2C^  century  US  Navy  destroyer,  through 
recent  ONR  Accelerated  Flydrodynamics  program.  Also,  simulations  for  naval  surface 
combatant  were  included  at  Gothenburg  2000  Workshop.  Levels  of  verification  and  validation 
and  overall  code  performance  demonstrated  that  CFDSFIIP-IOWA  was  among  the  best  codes  for 
the  combatant  test  case  (Wilson  et  al.,  2000). 

The  report  provides  documentation  of  the  modeling,  numerical  methods,  and  code 
architecture  for  CFDSHIP-IOWA  v3.03.  An  example  simulation  is  presented  to  demonstrate 
capabilities  and  concluding  remarks  are  provided.  This  report  and  a  suite  of  example  problems 
can  be  found  on  the  CFDSHIP-IOWA  website  http://www.iihr.uiowa.edu/~cfdship. 
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2.  CFD  PROCESS 

Overall  philosophy  for  the  CFD  process  is  given  in  this  section  as  a  set  of  procedures  to 
guide  engineers  and  scientists  through  the  process  of  modeling  fluid  flow  problems  using  a  CFD 
code.  Although  some  of  the  elements  of  the  CFD  process  are  relatively  straightforward, 
development  of  a  comprehensive  process  is  useful  for  training  non-expert  CFD  users, 
establishing  confidence  in  results  from  CFD  codes,  assessing  risks  in  the  use  of  CFD  results  in  a 
design  environment,  and  streamlining  the  task  of  obtaining  CFD  solutions  leading  to  reduced 
manpower  requirements.  As  described  in  the  following  paragraphs,  the  CFD  process  is 
composed  of  two  distinct  parts,  (/)  selection  or  development  of  a  general-purpose  CFD  code  and 
(ii)  use  of  the  CFD  code  for  solution  of  a  particular  flow  problem  of  interest.  In  general,  the 
former  occurs  only  at  infrequent  intervals  when  need  arises  to  make  large  shifts  in  technology, 
whereas  the  latter  must  be  followed  for  each  simulation. 

Development  of  any  general-purpose  CFD  code  has  several  common  elements. 
Specifically,  formulation  of  the  general  initial  and  boundary  value  problem  (IBVP)  which  is  to 
be  solved  numerically  using  a  CFD  code,  development  of  numerical  methods  for  approximate 
solution  of  the  IBVP,  and  documentation  of  the  CFD  code.  Key  issues  in  the  formulation  of  the 
IBVP  are  in  definition  of  the  scope  and  level  of  flow  description  (e.g.,  PANS,  LES,  DNS), 
selection  of  governing  partial  differential  equations  (PDEs)  and  physical  models  for  the  fluid 
flow,  and  selection  of  a  comprehensive  set  of  initial  and  boundary  conditions  required  to  solve  a 
wide  range  of  applications.  With  regard  to  numerical  methods,  key  issues  include  discretization 
of  the  continuous  modeled  PDEs,  initial  and  boundary  conditions,  development  of  numerical 
algorithms  for  solution  of  the  discretized  modeled  equations,  and  programming  and  testing  the 
algorithm  in  a  CED  code.  Einally,  documentation  of  the  CED  code  is  required  to  assist  users  in 
running  the  code.  Additional  documentation  may  be  required  to  assist  other  users  in  the 
development  of  new  models  or  numerical  methods. 

The  CED  process  for  simulation  of  a  fluid  flow  application  is  summarized  here  in  six 
distinctive  phases  as  shown  in  Eig.  1.  The  process  is  initiated  by  clearly  defining  the  purpose  for 
the  CED  simulation  and  the  acceptable  levels  of  CED  simulation  error  and  uncertainty  (e.g., 
prediction  of  vehicle  drag  with  validation  of  uncertainty  of  5%  over  a  specified  range  of 
Reynolds  numbers).  The  second  step  is  to  formulate  the  IBVP,  which  involves  definition  of  the 
required  governing  PDEs,  physical  models,  and  initial  and  boundary  conditions  for  the 
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application  of  interest.  In  addition,  the  flow  geometry,  domain,  and  coordinate  system  are 
defined  in  this  step.  As  an  aid  in  later  steps,  it  is  often  useful  to  eonstruct  sketehes,  sueh  as  those 
shown  in  Figs.  2  and  3,  whieh  summarize  geometry,  domain,  coordinate  system,  and  boundary 
eonditions.  It  is  assumed  that  the  seleeted  CFD  code  meets  the  requirements  for  the 


_ CFD  Process _ 

1.  Define  purpose  and  required  levels  for  V&V 

2.  Formulate  the  IB  VP 

•  Define  eontinuous  PDFs 

•  Define  physical  models 

•  Define  initial  and  boundary  eonditions 

•  Define  geometry,  domain,  coordinate  system 

3.  Plan  simulation  matrix 

4.  Create  input  files 

•  Generate  grid(s) 

•  Prescribe  boundary  conditions 

•  Prescribe  initial  conditions 

•  Select  flow  conditions 

•  Select  models 

•  Seleet  numerieal  parameters  and  post-proeessing  variables 

5.  Execute  CFD  code 

6.  Post-process  and  doeument  results 

Figure  1 .  CFD  Proeess 


particular  application  as  defined  in  this  step  or  ean  be  modified  to  do  so  with  an  aeeeptable  level 
of  effort.  If  further  code  development  is  required,  issues  with  souree  code  architecture  and 
availability  become  important  as  discussed  in  Section  4.  These  issues  may  impaet  whether  the 
user  selects  a  commercial  or  research  code. 

The  third  step  involves  planning  of  the  simulation  matrix  (i.e.,  number  of  simulations 
required  to  study  desired  range  of  flow  eonditions).  Depending  on  the  purpose  of  the  CFD 
simulation  and  the  environment,  some  or  all  cases  may  be  selected  to  estimate  simulation  error 
and  uneertainty  through  eomparison  of  the  simulation  results  with  available  benchmark  data. 
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Figure  2.  Open-Water  Propeller  Model  P5168; 
Computational  domain  and  boundary  conditions. 


Figure  3.  Surface  Combatant  Model  5415; 
Computational  domain  and  boundary  conditions. 


6 


Recommended  verification  and  validation  procedures  are  summarized  in  Section  7  to  accomplish 
this  task  and  involve  obtaining  solutions  with  multiple  grids  and  time  steps  with  systematic 
refinement.  As  a  result,  the  simulations  used  for  grid  and  time  step  studies  and  available 
benchmark  data  are  identified  in  this  step. 

The  fourth  step  of  the  CFD  process  involves  creation  of  the  CFD  code  input  files.  For 
CFDSFIIP-IOWA,  at  least  three  input  files  are  required  to  run  the  code  and  to  specify  the  grid 
coordinates  (FNAMEG),  boundary  conditions  (FNAMEO .  bcs),  and  runtime  parameters 
(cf  d_ship  .  nml).  The  purpose,  creation,  and  format  of  the  input  files  are  described  in  detail 
in  Section  6.  The  fifth  step  is  to  execute  the  CFD  code,  which  involves  selecting  a  computer 
platform  (e.g.,  serial  or  vector  machine),  compiling  the  source  code  to  produce  an  executable  if 
necessary,  and  running  the  CFD  code  as  described  in  Section  5.4. 

The  sixth  and  final  step  is  the  post-processing  and  documentation  of  results.  In  this  step, 
any  secondary  variables  of  interest  (e.g.,  vorticity,  turbulence  energy  budget)  are  computed  from 
the  primary  results  (i.e.,  velocity,  pressure,  turbulence  quantities).  If  the  results  are  unsteady, 
running  statistics  for  computation  of  variable  mean  and  variance  may  be  collected,  or  a  Fourier 
analysis  may  be  conducted  to  determine  the  harmonic  content  of  the  results.  Flow  visualization 
of  the  results  may  be  useful  in  interpreting  the  results  and  understanding  the  flow  phenomena.  If 
an  assessment  of  simulation  error  and  uncertainty  is  required,  a  verification  and  validation 
analysis  is  performed  at  this  phase  following  the  procedures  given  in  Section  7.  Finally,  the 
results  are  documented  in  an  appropriate  format. 

3.  MODELING 

As  mentioned  in  the  previous  section,  formulation  of  the  IBVP,  which  is  the  second  step 
of  the  CFD  process,  requires  definition  of  the  governing  equations,  physical  models,  and  initial 
and  boundary  conditions.  Flere,  this  definition  is  provided  through  a  detailed  presentation  of  the 
unsteady  RANS  equations,  turbulence  closure  and  model  equations,  initial  and  boundary 
conditions,  and  free-surface  and  body-force-propulsor  models. 
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3.1.  Governing  equations 


High-fidelity  resolution  of  a  certain  portion  of  the  frequency  spectrum  of  unsteady  flow 
physics  and  resultant  fluid  forces  and  moments  can  be  obtained  using  unsteady  RANS.  In  the 
context  of  a  triple  decomposition,  an  instantaneous  flow  quantity  can  be  written  as 

/(x,  t)  =  f{x)  +  f''{x,t)  +  f'{x,t)  =  F{x,t)  +  f'{x,t)  (1) 

where  /(x)  is  the  mean,  f"{x,t)  are  the  organized,  or  deterministic,  fluctuations,  and  f'{x,t) 
are  the  turbulent,  or  random,  fluctuations.  It  is  assumed  that  the  RANS  equations  solve  for 
F{x,t)  =  /(x)  +  f"{x,t)  and  that  the  Reynolds-averaging  process  is  based  upon  a  time  interval 
large  enough  to  average  out  f'{x,t)  but  also  small  enough  to  capture  f"{x,t)  .  This  implies  that 
the  frequencies  of  f"{x,t)  lay  sufficiently  outside  the  spectrum  of  turbulence  and  the  effect  of 
turbulence  upon  F(x,  t)  can  be  modeled  as  Reynolds  stresses.  This  also  implies  that  for  unsteady 
RANS,  all  turbulent  production  and  dissipation  is  subgrid  scale. 

The  code  has  been  formulated  to  solve  the  RANS  equations  in  either  Cartesian  or 
cylindrical-polar  base  coordinate  systems.  In  addition,  both  systems  may  be  in  either  absolute 
(i.e.,  earth-fixed  inertial)  or  relative  non-inertial  (i.e.,  fixed  to  a  moving  body)  reference  frames. 
Available  options  with  corresponding  input  parameter  values  are  listed  in  Table  1. 


Table  1.  Coordinate  system  options 


icoord 

Description 

Equations  solved 

1 

Cartesian,  absolute  frame 

(2)  -  (3) 

2 

Cartesian,  non-inertial  relative  frame 

(13) 

3 

Cylindrieal,  absolute  frame 

(4)  -  (8) 

4 

Cylindrieal,  non-inertial  relative  frame 

(4)  -  (7),  (14) 

For  Cartesian  coordinates  the  continuous  continuity  and  momentum  equations  in 
nondimensional  tensor  form  are 


^  =  0 

dUf  dU,  dp  1  d  - 

- -  +  Ui - ^  = - —  + - ^ - +  /, 

dt  dx.  dx,  Re  dxdx.  dx. 

J  ^  J  J  J 


(2) 

(3) 
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where  Ui  =  (U,  V,  W)  are  the  Reynolds-averaged  veloeity  eomponents,  x/  =  (x,  y,  z)  are  the 


independent  eoordinate  directions, 


P  = 


^  p-p, 

I  pul 


\ 

+  z/Fr" 

y 


is  the  piezometric  pressure 


coefficient,  u^Uj  are  the  Reynolds  stresses  which  are  a  two-point  correlation  of  the  turbulent 
fluctuations  ui,  /J  is  the  non-dimensional  body- force  vector  (^=  pUl^  where  is  a  force 

per  unit  volume  which  represents  the  effect  of  the  propeller,  Fr  =  U(^I the  Froude 
number,  and  Re  =  UoL!  v  is  the  Reynolds  number.  All  equations  are  nondimensionalized  by 


reference  velocity  Uq,  length  L,  and  density  p. 


For  cylindrical-polar  coordinates  the  continuity  and  momentum  equations  in 
nondimensional  vector  form  are 


dU  ^  \  d(rV) 
dx  r  dr 


1  dW 
r  dO 


=  0 


(4) 


DU 

Dt 


^  + J_ 

dx  Re 


d  —  d  —  15  — 1  1  * 

—  uu - uv - uw>-\ —  f 

dx  dr  rdP  \  p 


(5) 


DV  _  dp  ^  \  I  2j,  2  dW  V] 

Dt  r  dr  Re  \  r^  dO  r^j 


(6) 


DW  VW  \  dp  \  2  5F  W] 

Dt  r  r  dO  Re\  r^  dO  r^J 


(7) 


where  Ui  =  (U,  V,  W)  are  the  velocity  components  in  the  axial,  radial,  and  circumferential  (x,r,0) 
directions  and  the  operators  D/Dt  and  are  defined  as 


D  d  5  5  W  d 

Dt  dt  dx  dr  r  dO 

^2  5"  5'  1  5  1  5' 

V  = - 1 - 1 - 1 - 

dx^  dr^  r  dr  r^  dO^ 


(8) 


In  the  absolute  frame,  body  motions  are  resolved  by  time-dependent  grid  motions,  which, 
as  will  be  explained  in  section  3,  are  accounted  for  in  the  transformation  from  physical  to 
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computational  coordinates.  An  obvious  consequence  of  moving-body  simulations  in  the  absolute 
frame  is  that,  exeept  for  the  simple  cases  of  inertial  motions  such  as  steady  translation,  they  are 
inherently  unsteady. 

Relative  frame  formulations,  on  the  other  hand,  provide  capability  for  steady  simulation 
of  some  simple,  but  important,  non-inertial  motions.  Shown  in  Figure  4  are  two  examples: 
constant-speed  rotating  machinery  with  circumferentially  uniform  inflow,  such  as  a  propeller; 
and  a  ship  in  a  constant  radius  turn. 


(a)  Propeller  rotating  about  x-axis 


(b)  Ship  in  constant-radius-tum  about  z-axis 


Figure  4.  Example  relative-frame  applications. 


Transformation  to  relative  frame  is  straightforward.  The  acceleration  term  on  the  left  hand  side 
(LHS)  of  (3)  and  (5)-(7)  is  replaced  with  the  following  general  expression 


D\] 

Dt 


- 1 - ^-1 - xr  -l-nx(nxr  )-l-2nxU 

Dt  dt^  dt  ^  ’ 


(9) 


where  r  is  the  displacement  vector  of  the  relative  frame  with  respect  to  the  absolute  frame,  Q  = 
(rUjc,  o>y,  a^z)  is  the  angular  velocity  vector,  and  r'  and  U'  are  the  coordinate  and  velocity  vectors 
in  the  relative  frame,  respectively.  The  terms  on  the  right-hand-side  of  (9)  represent  non-inertial, 
linear,  tangential,  centripetal,  and  coriolis  accelerations  and  can  accommodate  general  6-degree- 
of- freedom  motions. 
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Currently,  however,  relative-frame  motion  in  Cartesian  eoordinates  is  restrieted  to  steady 
rotation  about  either  the  x-  or  z-axes.  For  these  simple  eases,  (9)  reduees  to  the  following 

-co^y'  -  2co^W'  -  (o]x:  +  2(oTJ' 

+  2(0/' 

where  (x',y',z')  and  {/' /' ,W'^  are  the  eoordinates  and  veloeity  eomponents  in  the  relative 
frame.  As  shown  in  Figure  4,  the  relationship  between  frames  is  a  funetion  of  time,  defined  by 
the  angle  P(t)  or  a(t).  In  eylindrieal  eoordinates,  motion  is  restrieted  to  steady  rotation  about  the 
x-axis,  whieh  reduees  (9)  to 

0 

-(oy-2(oJ¥' 

+2(0/' 

In  addition  to  modifying  the  aeeeleration  terms,  the  initial  and  boundary  eonditions  must  be 
transformed  into  relative  frame.  This  results  in  a  large  solid-body  rotation  of  the  free-stream 
veloeity  and  is  the  usual  approaeh  to  formulating  relative-frame  eodes  (e.g.,  Chen,  2000). 

An  alternative  approaeh  is  used  here  whieh  has  the  benefits  of  removing  the  solid-body 
rotation,  moving  most  of  the  non-inertial  terms  from  the  souree-term  on  the  RHS  to  the 
eonveetive  terms  on  the  LHS  of  (3)  and  (5)-(7),  and  simplifying  oaleulation  of  vortieity  and  wall- 
shear  stress,  sueh  that  the  same  algorithms  may  be  used  for  either  referenee  frame.  To  derive 
sueh  a  system  of  equations,  a  new  relative-frame  veloeity  veetor  {/" /" ,W"^  is  defined  with  the 


D\]  DV 

- — - h 

Dt  Dt 


DV  D\]' 

- — - h 

Dt  Dt 


solid-body  rotation  removed.  In  Cartesian  eoordinates 


-(O^Z  +  (O^X 


(12) 


Rearranging  (12)  for  (U'/',W'),  substituting  into  (2),  (3),  and  (10),  and  oolleeting  terms 
provides  RANS  equations  in  terms  of  {U"/",W") 
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dU"  dV"  dW"  ^ 

- + - + - =  0 

5x'  ej'  5z' 

du"  ,.du"  ,  ,.du"  ,.du" 

dt  ^  ’  dx'  ^  ’  dy'  ^  ’  dz' 

dp  1  ey  d  — 

= - ^ - MM,.  +  /.  -coY 

dx  Re  dxjdxj  dx'j 

dv"  ,yv"  ,  ,yv"  ,yv" 

- +  (U  -co^y) — r  +  (t/  -co^z  +co^x) — -  +  (U  — r  (13) 

dt  ^  ’  dx'  ^  ’  dy'  ^  ’  dz' 

dp  1  d^"  d  —  * 

= - —-I - VM  ,  +  fu  -  coW  +  coU 

dy'  Re  dx'dx'.  dx'. 

dw"  ,yw"  ,  ,yw"  ,yw" 

——  +  {U  -o}^y)—-i-  +  {U  -co^z  +o}^x)—-;-  +  {U  -a^y)—-;- 
dt  dx  dy  dz 

dp  1  d^W"  d  - 

=  — —  ^ - wu ,  +  fu  +  coy 

dz'  Re  dx'jdx'j  dx'j 

As  a  result,  all  of  the  eentripetal  and  half  the  Coriolis  terms  have  been  effeetively  moved  to  the 
LHS  of  (13)  in  the  form  of  modified  eonveetive  veloeities.  For  oylindrieal  eoordinates,  similar 
transformation  may  be  performed.  However,  in  this  ease,  all  of  the  non-inertial  terms  are  moved 
to  the  LHS  sueh  that  the  governing  equations  are  exaetly  the  same  as  (4)-(7)  with 
(Uy,JV)  =  (U"y",fV") ,  exeept  for  the  following  modifieation  of  the  substantial  derivative 

(14) 

Dt  dt  dx’  dr’  r'  d  O' 


It  is  noted  that  this  approaeh  is  essentially  the  same  as  the  Arbitrary  Lagrangian  Eulerian  (ALE) 
methods  (Hirt  et  al.,  1974)  wherein  the  eonveetive  veloeity  is  defined  as  the  differenee  between 
fluid  partiele  veloeity  and  the  mesh  veloeity. 

Although  there  are  numerous  subroutines  that  have  eoordinate-system  dependent  logie 
(e.g.,  boundary  eondition  formulation,  ealeulation  of  wall-proximity  funetions  and  wall-shear 
stress,  and  eoordinate  transformation  ineluding  grid-veloeity  terms),  users  only  are  required  to 
speeify  eonsistent  eoordinate  system  (icoord) ,  flow  eonditions  (agvx,  agvy,  agvz)  , 
and  boundary  eonditions,  the  latter  of  whieh  are  diseussed  in  Seetion  4.6. 
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3.2.  Turbulence 


CFDSHIP-IOWA  is  designed  to  use  a  linear  closure  model  where  the  Reynolds  stresses 
are  directly  related  to  the  mean  rate-of-strain  through  an  isotropic  eddy  viscosity  Vt.  In  Cartesian 
coordinates,  it  is  written  as 


-u,Uj  =  C 


dU, 


-  + 


dUj 

dX: 


--S^k 
3  ^ 


(15) 


where  dij  is  the  Kronecker  delta  and 


vv  +  ww) 


is  the  turbulent  kinetic  energy. 


In  cylindrical  coordinates, 


-MV  =  V, 


-uw  =  V. 


-vw  =  V, 


-MM  =  V, 


'dU 

y  dr  y 


--dyk 

3  ' 


--5yk 
3  " 


( 1  dU  dW'] 

- 1 - 

\r  d0  dx  J 

"idV  dW  W'^ 

^r  do  dr  r  j 


f 


--5± 


-vv  =  V, 


2^ 

V  dx , 


'^BV 

V  dr 


-ww  =  V, 


f2dW 


r  dO 


r  J 


(16) 


For  unsteady  flow,  equations  (15)  and  (16)  are  quasi-steady  relationships  where  it  is  assumed 
that  -u^Uj  responds  instantaneously  to  the  mean  rate  of  strain. 

Substituting  (15)  for  the  Reynolds-stress  term  in  (3),  the  momentum  equations  in 
Cartesian  coordinates  become 


dt  ^  dxj 


dP 

- h 

dx,  R, 


1  d^U,  dv, 


U,  dXjdXj 


-  +  - 


dx, 


dU,  dU, 


■  4"  - 
dx,  dx, 


+  fb, 


(17) 


where 


P  =  p  +  —k 
3 


(18) 
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(19) 


1  1 

- =  —  +  v, 

Ru  Re 


The  same  can  be  done  for  cylindrical  coordinates  where  (16)  is  substituted  into  (5), (6),  and  (7) 


Dt  dx  R, 


w, 

dv,  dU  dv, 
-  +  - 


+  2 

[  dx  dx 
DV 


dr 


( dU 
dr  dx 


15k 


+  — 


r  dO 


1  dU  dW 
- +  - 


Dt 


dP  1 

-2o)W-0)^r  = - + - W^V- 


r  dO  dx  j ^ 
2  dW  V 


1 

■+-4 

p 


dv. 


5x 

DW  VW 


(dU_  ^ 
dr  dx  y 


+  2 


Rv,  - 

dv,  dV  1  dv, 
-  +  — 


dr  dr  r  dO 


r  de  r\ 
dW 

r  80  dr  r  j 


1  y-.* 

■+-4 

p 


Dt 


.  T.  1 

—  +  2(oV  = - H - 

r  r  dO  R„ 


r"  de  r^ 


+  < 


dv. 


dx 


f  1  dU  dW^ 
r  de  dx 


+  - 


dv. 


dr 


dW 

^r  de  dr  r  j 


15k 


+  — 


r  de 


2dW  2,, 

- +  -V 

r  de  r  y 


1 

■+-4 

p 


(20) 


(21) 


(22) 


In  CFDSHIP-IOWA,  eddy  viscosity  can  be  calculated  using  one  of  two  models:  Baldwin- 
Lomax  or  the  blended  k-co/k-s  (BKW),  including  an  option  for  shear-stress  transport  (SST) 
model  (Menter,  1994).  The  turbulence  model  and  options  are  selected  using  the  input  parameters 
itm  and  itm_s witch  as  described  in  Section  6. 

For  the  Baldwin  Lomax  model  (itm=l),  details  have  previously  been  documented  in 
Stern  et  al.  (1996).  However,  as  illustrated  in  Paterson  and  Sinkovits  (1999),  care  must  be 
exercised  when  using  BL  for  geometries  with  multiple  no-slip  surfaces  and  multi-block  grid 
systems  since  the  turbulent  length  scale  is  calculated  as  a  weighted  average  based  upon  the  wall 
distance  to  each  no-slip  surface  in  a  given  block.  Since  a  search  process  through  all  blocks  is 
not  performed,  the  length  scale  may  be  incorrect  given  certain  blocking  topologies.  Therefore, 
BL  is  not  recommended  for  general  application  to  complex  geometries  and/or  grid  systems. 

The  BKW  model  (itm=2)  has  proven  to  be  robust,  applicable  to  complex  geometries  and 
flows,  and  fairly  accurate.  In  nearly  all  circumstances,  it  is  superior  to  A:-^  models,  which  require 
complicated  near-wall  models  that  are  difficult  to  implement  in  a  general  fashion.  The 
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governing  equations  for  the  eddy  viscosity  Vt,  turbulent  kinetic  energy  k,  and  the  turbulent 
specific  dissipation  rate  co  are  as  follows, 


(23) 


The  blending  function  Fj  was  designed  to  be  1  in  the  sublayer  and  logarithmic  regions  of 
boundary  layers  and  gradually  switch  to  zero  in  the  wake  region  to  take  advantage  of  the 
strengths  of  the  k-co  and  k-f  models,  i.e.,  k-co  does  not  require  near-wall  damping  functions  and 
uses  simple  Dirichlet  boundary  conditions  and  the  k-s  does  not  exhibit  sensitivity  to  the  level  of 
free-stream  turbulence  as  does  the  k-co  model.  The  distance  to  the  nearest  no-slip  surface  d  is 
required  for  calculation  of  Fi  and  the  model  constants  are  calculated  locally  as  a  weighted 
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average,  i.e.,  </)  =  F^(/)^+{l- F^^(/)^  where  are  the  standard  k-co  model  and  (p2  are  the 


transformed  A:-^  model  constants  in  Table  2. 

Table  2.  Blended  k-oj/k-s  model  constants. 


4> 

(F 

(j)2 

SST 

Ok 

0.5 

1.0 

0.85 

^co 

0.5 

0.856 

0.5 

P 

0.075 

0.0828 

0.075 

P* 

0.09 

0.09 

0.09 

K 

0.41 

0.41 

0.41 

Y 

0.0553 

0.04403 

0.0553 

In  addition  to  the  standard  BKW  model,  a  SST  model  (M enter,  1994)  is  included  as  a  user 
specified  option.  The  SST  model  accounts  for  transport  of  the  principal  turbulent  stresses  and 
has  shown  improved  results  for  flows  with  adverse  pressure  gradients.  The  SST  model  is 
identical  to  the  standard  model  except  for  a  change  in  <Jk  as  shown  in  Table  2  and  the  definition 
of  eddy  viscosity 


V,  = 


0.3  lA: 


max  (0.3 10,0^2) 


Fj  =  tanh 


max 


lyfk  SOOu 


0.09a>y  y^co 


(29) 


where  f2  is  the  absolute  value  of  the  vorticity,  and  y  is  the  distance  to  the  nearest  wall.  This 
effectively  reduces  eddy  viscosity  in  regions  of  high  off-body  vorticity  such  as  that  found  in 
separated  flow  or  in  a  tip  vortex. 


3.3.  Initial  and  Boundary  Conditions 

Formulation  of  an  IB  VP  requires  mathematical  derivation  of  the  initial  and  boundary 
conditions  for  each  dependent  variable  and  for  each  type  of  condition  that  is  to  be  simulated. 
Complete  presentation  of  the  available  palette  of  conditions  and  underlying  numerics  in 
CFDSHIP-IOWA  is  deferred  until  Sections  4.5  and  4.6.  Here,  a  brief  discussion  of  the  initial 
and  boundary  condition  modeling  is  provided. 
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Since  CFDSHIP-IOWA  can  be  run  in  two  different  modes,  i.e.,  steady  (time-asymptotic) 
and  unsteady  (time-accurate),  the  initial  eonditions  serve  two  purposes.  For  steady-flow 
simulations,  the  initial  conditions  provide  the  zeroth  iteration  to  an  iterative  scheme  and  can  be 
fairly  crude.  Usually,  out  of  eonvenience,  free-stream  conditions  are  used  throughout  the 
domain.  For  unsteady  flow,  on  the  other  hand,  the  initial  conditions  serve  as  the  solution  at 
time=0.0  and  should  therefore  satisfy  the  governing  equations  at  this  time.  General  speeification 
for  arbitrary  geometries  and  eonditions  is  nearly  impossible;  therefore,  the  available  initial 
condition  corresponds  to  a  start  from  rest  (i.e.,  U  =  V  =  W  =  P  =  0.0,  k  =  k/st,  co  =  co/st)  with  a 
cubic  polynomial  acceleration  to  ship  speed.  Details  of  numerical  implementation  are  provided 
in  Section  4.5. 

As  with  most  CFD  codes,  there  are  numerous  BC  types  which,  for  convenience,  can  be 
grouped  into  domain  truncation  boundaries,  physieal  boundaries,  and  computational  boundaries. 
Physieal  BC’s  are  due  to  solid  surfaces  or  water-air  free  surfaee,  the  latter  of  which  is  described 
in  the  next  section.  For  external-flow  hydrodynamics,  an  infinite  unbounded  fluid  often 
represents  the  physical  domain.  This  requires  that  the  computational  domain  be  truneated  to  a 
size  that  ean  be  eeonomically  filled  with  grid  points,  but  has  no  influenee  on  the  computed 
solution.  Computational  boundaries  are  due  to  grid  topology,  modeling  assumptions,  and  multi- 
bloek  domain  decomposition.  Available  options  in  CFDSHIP-IOWA  are  listed  in  Table  7  and 
will  be  discussed  in  more  complete  detail  in  Section  4.6. 


3.4.  Free-surface 

CFDSHIP-IOWA  uses  a  surface  traeking  approach  for  modeling  the  free  surface.  The 
kinematic  free-surface  boundary  condition  (KFSBC)  is  used  to  eompute  the  evolution  of  the  free 
surface,  while  the  dynamic  free-surface  boundary  condition  (DFSBC)  provides  boundary 
eonditions  for  veloeity  and  pressure.  Considering  the  KFSBC,  the  requirement  that  the  wave 
elevation  ^"be  a  stream  surfaee  is  satisfied  by  the  condition 

D(^-z) 


Dt 


=  0 


Expanding  equation  (30)  gives  a  continuous  2D  hyperbolic  PDE  for 

^  +  U^  +  V^-W  =  0 
Ft  dx  dy 


(30) 


(31) 
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At  the  intersection  of  free-  and  no-slip  surfaces  (i.e.,  the  contact  line),  equation  (31)  becomes 
singular  when  the  contact  line  is  in  motion  but  the  fluid  velocity  is  zero  due  to  the  viscous  no-slip 
boundary  condition.  This  problem  is  overcome  through  the  use  of  an  approximate  contact  line 
model  where  a  small  near-wall  region  is  “blanked  out”  when  solving  equation  (31)  and  the 
solution  in  this  region  is  linearly  extrapolated  from  the  interior  of  the  domain.  The  numerical 
method  for  solution  of  the  KFSBC  given  by  equation  (31)  is  presented  in  Section  4.4. 

The  DFSBC  requires  that  the  normal  and  tangential  stresses  are  continuous  across  the 
free-surface 


(32) 


where  nj  is  the  unit  normal 
=  -pS.j  +  Re^'  +  dU j  jdx^  ^  -  upj 


vector  to  the  free  surface  and  rij 
and  Tij*  are  the  fluid-  and  external-stress  tensors, 


respectively.  Although  not  included  in  this  presentation,  the  effects  of  air  and  surface  tension 
can  be  included  through  the  external-stress  tensor  rij* .  The  following  approximations  were  used 
to  obtain  free-surface  boundary  conditions  from  equation  (32):  (/)  the  external  stress  and  surface 
tension  are  assumed  zero  so  that  T^rij  =  0 ;  and  (zz)  the  gradients  of  the  free  surface  and  normal 


velocity  in  the  tangential  directions  are  assumed  small  (i.e.,  d(^ jdx  0  0,5^ jdy  □  0 ,  dWjdx  □  0 , 
and  dWjdyU  0).  Under  these  assumptions,  expansion  of  equation  (32)  gives  the  following 
approximate  dynamic  boundary  conditions  for  pressure  and  tangential  velocity  components 

<«) 

Lastly,  a  zero-gradient  condition  is  used  for  W,  which  is  consistent  with  the  approximations 
employed  for  the  dynamic  condition 


dW 

dz 


=  0 


(35) 


5.5.  Body-force  propulsor 


The  momentum  equations  (3)  include  a  body-force  term  ,  which  may  be  used  to  model 
the  effects  of  a  propulsor  without  resolving  the  detailed  blade  flow.  There  are  numerous 
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approaches  to  calculating  including  simple  prescribed  distributions,  which  recover  the  total 

thrust  and  torque,  to  more  sophisticated  methods  which  use  a  propeller  performance  code  in  an 
interactive  fashion  with  the  RANS  solver  to  capture  propeller-hull  interaction  and  to  distribute 
according  to  the  actual  blade  loading  (e.g..  Stem  et  al.,  1994).  For  the  latter,  custom  interface 


software  must  be  developed  to  extract  effective  wake  from  RANS  solution  and  to  produce  /* 


calculated  by  performance  code.  This  is  not  provided  with  CFDSHIP-IOWA,  but  can  be  easily 
developed  by  experienced  users  with  access  to  a  propeller  performance  code. 

CFDSHIP-IOWA  does,  however,  include  a  prescribed  axisymmetric  body  force  with 
axial  and  tangential  components  (Stem  et  al.,  1988).  The  radial  distribution  of  forces  is  based 
the  Hough  and  Ordway  circulation  distribution  (Hough  and  Ordway,  1964)  which  has  zero 
loading  at  the  root  and  tip.  Therein, 


fbx  =  A/ yj\-r* 
f  -A 

Jb9  ~  ^9 


(\-R„y+R, 


(36) 


where 


,  r/RP-RH 

r  =  — - 

1-RH 

r  =  ^(_y  -  Y_PR0P_CENTER)^  +[z-  Z_PR0P_CENTER)^ 

^  _ 105  (37) 

""  (dXPROP)  16(4-t3RH)(l-RH) 

^  Kq _ 105 

^  (dXPROP)J^  ;r(4-t3RH)(l-RH) 


and  where  Ct  and  Kq  are  the  thrust  and  torque  coefficients,  J  is  the  advance  coefficient,  RP  is  the 
propeller  radius  non-dimensionalized  by  ship  length,  RH  is  the  hub  radius  in  decimal  percent  of 
RP,  and  DXPROP  is  the  mean  chord  length  projected  into  the  x-z  plane.  As  derived,  these  forces 
are  defined  over  an  "actuator  cylinder"  with  volume  defined  by  RP,  RH,  and  DXPROP,  i.e., 

;r^RP^  (l-RH^jj  DXPROP .  Integration  of  the  body  forces  (36)  over  this  analytical  volume 
exactly  recovers  the  prescribed  thmst  and  torque. 
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T  =  plJU 0 1 / 1/ 1  ”  fbjdxdOdr 
Q  =  pL^Ul  j  ^  f  j  "  dxd0 dr 

JRh  jo  Jx^ 


(38) 


Since  curvilinear  non-orthogonal  multi-block  grids  are  typically  used,  implementation  of 
(36)  requires  several  issues  to  be  addressed.  A  vertex-based  search  algorithm  is  used  to 
determine  which  grid-point  control  volumes  are  within  the  actuator  cylinder. Figure  5  shows  an 
example  for  commercial  ship  geometry,  Esso  Osaka.  In  this  simulation,  an  0-grid  topology  is 
used  to  wrap  around  the  stem.  The  search  algorithm  identifies  814  total  cells  in  two  blocks  that 
lie  within  the  cylinder  and  upon  integration  give  an  approximate  volume  of  the  cylinder,  i.e., 
prescribed  volume  =  1.6979x10'^  and  integrated  volume  =  1.6714x10'^.  Given  this  1.6%  error  in 
volume,  total  thrust  and  torque  in  (38)  is  not  recovered.  Therefore,  magnitude  of  body  forces  in 
(36)  are  uniformly  scaled  by  the  volume  error  such  that  the  integrated  total  force  is  equal  to  that 
which  is  prescribed.  Finally,  it  should  be  noted  that  all  body- force  algorithms  are  designed  to 
only  work  with  Cartesian  coordinates  (icoord=l  or  2). 


4.  NUMERICAL  METHODS 

This  section  describes  the  numerical  methods  used  in  CFDSHIP-IOWA  and  includes  a 
discussion  of  the  coordinate  transformation,  discretization  scheme,  free-surface  solver  and 
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adaptive  gridding,  RANS  solution  algorithm  and  pressure-Poisson  equation,  initial  and  boundary 
eonditions,  Chimera  overset  gridding,  and  ealeulation  of  forces  and  moments. 


4.1.  Coordinate  transformation 

The  continuous  governing  equations  are  transformed  from  the  physical  domain  in  either 
Cartesian  (x,y,z,t)  or  cylindrical-polar  {x,r,6,t)  coordinates  into  the  computational  domain  in  non- 
orthogonal  curvilinear  coordinates  {^,  rj,  f,  r).  A  partial  transformation  is  used  in  which  only  the 
independent  variables  are  transformed,  leaving  the  velocity  components  Ui  in  the  base 
coordinates.  The  transformation  relations  are 


J  ' 


vV  = 


1  5 


J  df 


Jg 


ij  ^(t> 


=  g 


d(j)  _d(j)  1  dx.  dcj) 

dt  dr  J  ‘  dt 


(39) 

(40) 

(41) 

(42) 


where  qt  represents  the  components  of  an  arbitrary  vector  q{xi).  The  geometric  coefficients  bj 


and  ,  the  Jacobian  J ,  and  /'  are  functions  of  coordinates  only  and  are  defined  for  Cartesian 
grids  as 


(43) 

(44) 


J  = 


x^  x^  x^ 
Tf  T,  Tf 


(45) 


(46) 
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where  is  the  permutation  symbol  with  Imn  eyelie.  The  grid-veloeity  terms  — -  in  (42), 

dt 

which  are  used  only  for  unsteady  flows  in  an  absolute  frame  of  reference,  are  calculated  directly 
using  finite  difference  expressions,  as  given  in  the  following  section. 

Using  the  transformation  relations,  the  continuity  (2)  and  momentum  (17)  equations  are 
written  as 


1  5 


(*/£/,)  =  0 


dr 


—  + a. 


dU..  I  rk  dP  1 


0^' 


-L  =  -  —  b'" 


■  +  ■ 


j  ‘  R 


-g 


^km 


-eff 


d^U. 

0^^ 


■  + 


'-4‘ 


J  "  0^' 


1  dU. 

J^m  ] 


J  '  0^” 


+  /w 


where 


k  ^  ul 
an  =—b, 


1,™0U  0X,^ 


J  '  J  '  0^"  0r ) 


R 


eff 


(47) 

(48) 


(49) 


Note  that  the  convective-term  coefficient  in  (49)  contains  contributions  from  both  the  linear 


Reynolds-stress  closure  (15)  and  the  grid  velocity  (42),  the  latter  of  which  introduces  non- inertial 
accelerations  due  to  body  motions. 

Splitting  the  viscous  term  into  normal  and  cross  components  and  rearranging  gives  the 
continuous  form  of  the  momentum  equations  in  the  computational  domain 


dU,  ,  dU,  1  ,,  0'U,  1  - ,  dP 

- ^  — i - S  — ^ 

dr  0^'  R^ff  d^'d^‘  J  ‘  0^' 


(50) 


R.. 


-12 


g 


0'U, 

07W 


■+g 


13 


0'U, 

8^' 


-  +  g 
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0'U. 

0^^ 


+ 

( 1 

y 

yj  0rj 

+  fM  (51) 


4.2.  Discretization  scheme 

For  temporal  discretization  of  equations  for  k-co  (24),  KFSBC  (31),  and  momentum  (50), 
a  general  formula  for  an  Euler  backward  difference  is  given  by 

^  =  (52) 

0r  Ar 
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where  the  weights,  wtn,  wtm,  wtmm,  determine  the  order  of  the  differenee  expression  and  are  given 
in  Table  3  for  the  first-  and  seeond-order  formulations.  For  steady  state  and  time-aceurate 
solutions,  first-order  formulation  and  second-order  formulations  are  used,  respectively,  and  are 


Table  3.  Finite-difference  weights  for  temporal  discretization. 


Scheme 

itemp_order  wtn 

wt„ 

f‘  order 

1  1 

-1 

0 

2"“^  order 

2  3/2 

-2 

1/2 

specified  using  the  input  variable  itemp_order.  Eq.  (52)  is  also  used  to  compute  grid 

3x- 

velocity  terms  — -  in  Eq.  (42)  where  the  general  variable  (j)  is  replaced  with  the  grid  coordinates 
dt 

Xi. 

Eor  spatial  discretization  of  equations  for  k-co  (24),  KESBC  (31),  geometric  coefficients 
(43),  Jacobian  (45),  and  momentum  (50),  the  convective  (or  first  derivative)  terms  are  discretized 
with  the  following  higher-order  upwind  formula 

^  =  -\u‘\)Sl^  (53) 

where 

<!>  =  +  wj,  +  W (54) 

^  (55) 

Six  convective  schemes  are  available  in  CEDSHIP-IOWA  and  their  weighting  coefficients  are 
supplied  in  Table  4. 

Table  4.  Finite-difference  weights  for  spatial  discretization  of  convective  terms. 


Scheme 

ispat_order 

^mm 

w„ 

Wp 

Wpp 

1  order  upwind 

1 

0 

-1 

1 

0 

0 

2"“*  order  eentral 

2 

0 

-1/2 

0 

1/2 

0 

2"“*  order  upwind 

3 

1/2 

-2 

3/2 

0 

0 

2"“^  order  upwind  biased  (Quiek) 

4 

1/8 

-7/8 

3/8 

3/8 

0 

3'^^  order  upwind  biased 

5 

1/6 

-1 

1/2 

1/3 

0 

4*  order  eentral 

6 

1/12 

-2/3 

0 

2/3 

1/12 
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The  user  may  specify  the  order  of  accuracy  for  the  momentum  equations  and  the  KFSBC  using 
the  input  variable  ispat_order,  however,  2"‘*-order  upwind  is  sufficient  for  most  simulations. 
For  the  BKW,  equation  (24)  is  discretized,  by  default,  using  the  same  scheme  as  the  momentum 
equations.  However,  order-of-accuracy  can  be  set  to  a  lower-order  scheme  using  namelist 
variable  itm_spat_order.  To  maintain  stability,  it  is  occasionally  necessary  to  set  the 
turbulence  model  discretizaiton  to  H*-order  upwind.  Evaluation  of  the  transformation  relations 
in  equations  (43)  and  (45)  is  accomplished  using  the  2"‘^-order  central  scheme. 

The  viscous  terms  are  written  as 
(h 

TT—  =  G)^n,m  +  ®2„  +  (Ol  +  (Ol  ^,^2  (56) 

oq  og 

and  the  coefficients  are  supplied  in  Table  5.  Although  written  in  a  general  fashion,  the  viscous 
terms  are  preset  to  2”‘*-order  central,  and  selection  of  4*’'-order  scheme  is  not  accessible  by  input 
parameter. 

Table  5.  Finite-difference  weights  for  spatial  discretization  of  viscous  terms. 


Scheme 

Order 

CO^m 

a)2n 

a2p 

a2pp 

2”‘^  order  central 

0 

1 

-2 

1 

0 

4*  order  central 

4th 

-1/12 

16/12 

-30/12 

16/12 

-1/12 

Applying  the  temporal  and  spatial  discretizations  given  by  equations  (52),  (53),  and  (56) 
to  the  continuous  momentum  equations  (50)  gives  the  discrete  form  of  the  momentum  equations 


nb 


(57) 


where  and  Ant  denote  the  central  and  neighboring  coefficients  of  the  discretized  momentum 
equations,  respectively.  The  source  term  contains  velocities  from  the  previous  two  time 


steps  (n-1)  and  (n-2)  and  the  mixed  derivative  terms  (51),  the  latter  of  which  are  lagged  to  the 
previous  iteration. 


1 


Su,  =%,-  -  —  ) 

'  Ar 


(58) 
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4.3.  RANS  solution  algorithm  &  pressure  Poisson  equation 


The  pressure-implicit  split-operator  (PISO)  algorithm  for  solving  the  incompressible 
Navier-Stokes  equations  (Issa,  1985)  uses  a  predictor-corrector  approach  to  advance  the 
momentum  equation  while  enforcing  the  continuity  equation.  In  the  predictor  step,  the 
momentum  equation  (57)  is  advanced  implicitly  using  the  pressure  field  from  the  previous  time 
step 


nb 


dpn-^ 


(59) 


where  superscript  is  used  to  denote  advancement  to  an  intermediate  time  level. 
In  the  corrector  step,  the  velocity  is  updated  explicitly 


JAj,  Qe 


(60) 


using  a  pressure  obtained  from  a  derived  Poisson  equation  and  where  the  psuedo-velocity  is 
defined  as 


U.  =■ 


A. 


ijk 


s,-'Ea„,u; 


nb 


nb 


A  pressure -Poisson  equation  is  derived  by  taking  the  divergence  of  equation  (60) 


15  ■»»  Id  ■  Id 

-—b!UT=-—b!U,  --—b! 


j  d^^ 


j  d^^ 


j  d^^ 


1 


dP 


and  by  realizing  that  the  LHS  of  equation  (62)  goes  to  zero  upon  convergence 


1  5 


J  d^^ 


dP 


LA. 

J  d^^ 


biU, 


(61) 


(62) 


(63) 


Because  a  regular,  or  collocated,  grid  approach  is  used,  solution  of  equation  (63)  requires  special 
treatment  to  avoid  odd-even  decoupling.  Fourth-order  artificial  dissipation  is  implicitly  added  by 
taking  a  linear  combination  of  full-  and  half-cell  operators  (Sotiropoulos  and  Abdallah,  1992) 

{\-y)LP*  +  yLP*  +  NP*  (64) 


where  L  is  the  full-cell  formulation,/,  is  the  half-cell  formulation,  and  N  is  the  operator 
containing  mixed-derivative  terms 
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^=7^1  (“‘‘<5s)+4(“%)+4(‘'”4)}  (*5) 

i = 7{4  (“"4 ) + 4  (“”4 ) + 4  (“”4 ))  («) 

Af  =  i  {<?,  (a"4  +  a"4 )  +  4  ^  (<,='<?.  +  c,’% )}  (67) 

and  where 

<5..  !«  =  7(4, (68) 

4^  =  (^w/2 -^,-1/2)  (69) 


The  weighting  funetion  y  ranges  from  1  (i.e.,  most  dissipation  and  smooth  solutions)  to  0  (i.e.,  no 
dissipation,  but  prone  to  decoupling).  This  parameter  is  set  by  the  namelist  variable  gama_pr 
which  has  a  default  value  of  y  =  1.0.  Use  of  the  half-cell  operator  introduces  metrics  at  half-cell 
locations  which  are  computed  from  an  average  of  the  nodal  values.  Note  that  the  half-cell 
formulation  achieved  with  y=  1.0  is  essentially  the  same  as  the  Rhie  and  Chow  (1983) 
interpolation  method  for  avoiding  odd-even  decoupling. 


4. 4.  Free-surface  solver  and  adaptive  gridding 

The  tracking  approach  used  in  CFDSHIP-IOWA  for  modeling  the  free  surface  was 
presented  in  Section  3.3.  Therein,  the  DFSBC  was  used  to  provide  boundary  conditions  for 
velocity  and  pressure  as  given  by  equations  (33)-(35),  which  are  relatively  straightforward  to 
numerically  implement.  The  KFSBC  was  developed  by  requiring  that  the  free  surface  be  a 
stream  surface  resulting  in  a  2D  PDF  for  the  evolution  of  the  wave  elevation  ^ ,  as  given  by 
equation  (31).  The  numerical  method  for  solution  of  this  equation  will  be  presented  in  this 
section  and  closely  follows  the  approach  for  solution  of  the  RANS  equations  presented  in 
Sections  4.1  and  4.2. 

Equation  (31)  is  transformed  from  the  physical  domain  in  Cartesian  coordinates  (x,y,t) 
into  the  computational  domain  in  non-orthogonal  curvilinear  coordinates  using  a  reduced 

3D  form  (i.e.,  two  spatial  and  one  temporal  coordinate)  of  the  general  4D  transformation 
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presented  in  Seetion  4.1.  Using  the  transformation  relations,  the  eontinuous  KFSBC  in 
eomputational  spaee  is  given  by 

^  +  aJ-^-lF  =  0  (71) 

5r 

where  the  superseript  ‘A:’  is  summed  for  A=l,2  in  equation  (71)  and 


,  1  dx, 

a'i=-b':  U, - ^ 


The  temporal  term  in  equation  (72)  is  diseretized  using  an  Euler  baekward  differenee  as  given  by 
equation  (52),  while  the  eonveetive  term  is  diseretized  using  the  same  higher-order  upwind 
differenee  as  for  the  eonveetive  term  of  the  RANS  equations  and  is  given  by  equation  (53). 
Applying  the  temporal  and  spatial  diseretizations  to  the  eontinuous  KFSBC  (71)  gives 

=  (73) 

nb 

where  4^  and  [ehanged  ‘ij’  subseript  to  ‘nb’]  denotes  the  eentral  and  neighboring 

eoeffieients  of  the  diseretized  KFSBC.  The  souree  term  Si  eontains  the  vertieal  veloeity 
eomponent  W  and  the  wave  elevation  at  previous  time  steps. 

A  straightforward  solution  of  equation  (73)  often  leads  to  stability  problems  due  to 
several  faetors.  First,  as  diseussed  in  Seetion  2.2,  equation  (73)  is  singular  at  the  eontaet  line. 
Seeondly,  a  eombination  of  highly-elustered  near-wall  spaeing  required  for  turbulenee  models 
(i.e.,  on  the  order  of  10'^),  high-aspeet  ratio  grid  eells  (i.e.,  on  the  order  of  10^),  and  laek  of  either 
physieal  or  numerieal  dissipation  results  in  an  unstable  numerieal  system.  As  diseussed  in 
Seetion  2.2,  the  eontaet-line  is  modeled  by  “blanking  out”  the  solution  in  the  near  wall  region 
and  extrapolating  the  solution  from  the  interior.  The  blanking  distanee  is  set  by  the  variable 
wavblank  in  the  $FREE_SURFACE  namelist.  In  addition,  the  solution  is  filtered  with  a 
variable-order  high-bypass  filter  after  each  iteration 

Ci  =  +  2  (  ^,+1  +  Ci-\  )  +  2  (  C+2  +  Ci-2  )  ”'"  ^  (  ^(+3  +  '^1-3  )  ‘*) 

where  the  coefficients  a,  b,  c,  and  d  are  given  in  Table  6  for  second-,  fourth-,  and  sixth-order 
filters. 
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Table  6.  Filter  coefficients. 


Order 

a 

b 

c 

d 

V2 

V2 

0 

0 

qth 

5/8 

V2 

-1/8 

0 

11/16 

15/32 

-3/16 

1/32 

By  transforming  the  filter  from  physical  to  wave  number  space  (Lele,  1992),  it  can  be  shown  that 
the  4*  and  6*^  order  filters  remove  energy  at  the  highest  wave  number  ( k^sx  U  tt)  while  leaving 
the  low  wave  numbers  unchanged.  Since  the  order  fdter  is  overly  dissipative,  its  use  should 
be  avoided.  The  filter  coefficients  are  specified  in  the  boundary  condition  file  (if  sf  ilter=  2, 
4,  or  6)  and  are  default  to  the  4*-order  filter. 

The  discrete  KFSBC  given  by  equation  (73)  is  solved  on  all  block  faces  identified  as  a 
free-surface  boundary  using  the  iterative  solvers  as  described  in  Section  3.2.  After  the  filtering 
operation,  the  solution  is  used  to  conform  the  volume  grid  to  the  new  wave  elevation  through  the 
use  of  cubic-spline  interpolation  in  the  curvilinear  coordinate  of  the  original  grid  system.  This 
is  equivalent  to  a  “rigid-wire”  approach  where  the  grid  points  slide  along  the  ^  coordinate  as 
shown  in  Figure  6. 


Figure  6.  Dynamic  Free-Surface  Adaptive  Grid. 


In  general,  this  approach  is  fairly  robust,  but  is  susceptible  to  poor  grid  quality  if  there  is  large 
difference  between  the  original  and  conformed  grid  or  if  there  is  severe  geometry  changes  in  the 
,  or  girthwise,  direction. 
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4.5.  Initial  conditions 


Solution  of  the  IBVP  requires  initial  eonditions.  For  steady-flow  simulations,  the  time- 
marehing  proeess  serves  as  an  iteration  loop  and  the  initial  eonditions  provide  the  zeroth 
iteration.  Usually,  a  fairly  “poor”  initial  guess  (e.g.,  free-stream  at  all  grid  points,  exeept  for  no¬ 
slip  boundaries)  is  suffieient  sueh  that  the  algorithm  is  eapable  of  damping  out  any  initial 
transients.  For  unsteady  flow,  on  the  other  hand,  the  initial  eonditions  serve  as  the  solution  at 
time=0  and  initial  eonditions  whieh  do  not  satisfy  governing  equations  ean  prevent  the 
simulation  from  eonverging.  Exeluding  eustom  and/or  novel  unsteady  problems,  most 
applieations  ean  be  served  by  the  two  initial  eondition  options  available  in  the  eode. 

The  first  option  (mode=0)  sets  all  variables  to  uniform  free  stream,  i.e.,  the  dependent 
variables  have  the  following  values  f/=UINF,  F=VINF,  1U=WINF,  p=t),  k=kfst=lO'^,  a)=ci)fst=9.0, 

o 

and  V//i(=l.lxlO'  where  (UINF,  VINF,  WINF)  permits  speeifieation  of  free-stream  unit 


veetor.  For  steady  flow,  an  impulsive  start  is  used  where  veloeity  is  set  to  no-slip  values  at  the 
first  time  step  (iteration).  In  eontrast,  for  unsteady  flow,  the  no-slip  boundaries  are  smoothly 
ramped  from  free-stream  values  to  no-slip  values  using  a  eubie  polynomial.  For  time  < 
t  ime_ramp_end,  the  no-slip  boundaries  are  set  to  the  following 


U  =  U  iffp[\-ramp) 

V  =  Vjj^p{l-ramp) 
W  =  Wjj^p{\-ramp) 


ramp  =  -2 


time_ramp_end 


time_ramp_end 


(75) 


This  represents  a  smooth  aeeeleration  of  the  ship  from  rest  and  provides  initial  eonditions  that 
satisfy  continuity. 

The  second  option  (mode=l)  sets  initial  conditions  by  reading  a  restart  file.  While  a 
previous  simulation  typically  generates  this  Fortran  binary  file,  it  can  also  be  created  by  the  user 
to  prescribe  either  initial  conditions  or  boundary  conditions.  For  the  latter,  the  restart  file  must 
be  used  in  conjunction  with  boundary  condition  ibtyp  =14  which  is  described  below.  Format 
of  the  restart  file  is  described  in  Appendix  A.  6. 
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4.6  Boundary  conditions 


As  discussed  in  Section  2,  the  CFD  Process  requires  formulation  of  an  IBVP  where 
boundary  conditions  (BC)  must  be  specified  on  all  faces  of  the  computational  domain.  As  with 
most  CFD  codes,  there  are  numerous  BC  types,  which  for  discussion  here,  can  be  grouped  into 
domain  truncation  boundaries,  physical  boundaries,  and  eomputational  boundaries.  The 
formulation  of  each  BC  type  is  described  in  detail  and  guidanee  provided  on  when  and  how  eaeh 
BC  used. 

Twenty-six  different  BC  types  are  available  in  CFDSHIP-IOWA  and  are  summarized  in 
Table  7.  Eaeh  faee  of  each  mesh  must  be  speeified  and  ean  be  broken  into  an  arbitrary  number 
of  rectangular  patches  over  which  a  different  boundary  condition  can  be  applied.  Numerically, 
the  26  different  eonditions  eonsist  of  combinations  of  Dirichlet  and  Neumann  boundary 
conditions  for  the  nondimensional  flow  variables  U,  V,  W,  p ,  k,  a>  and  h.  Note  that  BC  for  Vt 
are  implemented  so  that  eddy-viseosity  gradients  in  equation  (17)  can  be  evaluated  near 
boundaries  without  changing  finite-difference  stencil.  Dirichlet  conditions  are  prescribed  values 
or  lagged  data  from  donor  regions  (i.e.,  periodie  or  multiblock).  Neumann  eonditions  are 
prescribed  gradients,  which  are  evaluated  using  one-sided  finite  differenees.  For  zero-gradient 
conditions,  two  functions  are  used  in  CFDSHIP-IOWA  to  evaluate  first  and  second  derivatives 


zero_fd  =  a-l-ibcord^(a-h) 


zero  sd  =  2a-h 


(76) 

(77) 


where  a  and  b  are  the  values  one  and  two  grid  points  inside  the  boundary  in  the  ibdir 
eoordinate  direction,  respectively.  Details  unique  to  each  BC  type  are  now  described. 
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Table  7.  Boundary  Conditions 


IBTYP 

Description 

u 

V 

w 

p 

k 

CO 

u 

c/5 

'C 

10 

Inlet 

UINF 

VINF 

WINF 

dPjd^,  =  0 

kfs,=lxl0-' 

K)fst=9.0 

Vt.fst 

-a 

p 

p 

o 

PQ 

11 

Exit 

=  0 

ev/ei/  =  0 

d^wld4f=0 

dP/d^,  =  0 

ayt/ai,  =  0 

da>/d^,  =  0 

dvjd^,=0 

1 

"cd 

O 

12 

Far- field  #1 

UINF 

dVjd^,  =0 

dWjd^,  =0 

0 

a/t/a#,  =  0 

a®/a|,  =  0 

ar./ai,.  =  0 

1 

H 

13 

Far-field#2 

UINF 

VINF 

WINF 

dPjd^,  =  0 

ai/ai,  =  0 

ary/ai,  =  0 

a^/ai,  =  0 

§ 

o 

Q 

14 

Prescribed 

* 

* 

* 

* 

* 

* 

* 

20 

Absolute-frame  no-slip 

0 

0 

0 

dPjd^,  =  0 

0 

6o/  Re  PAy‘ 

0 

c/5 

22 

Relative-frame  no-slip 

X 

i” 

z 

dP/d^,  =  0 

0 

6o/  Re  PAy‘ 

0 

T3 

§ 

27 

Impermeable  slip  (calculate 

Eq.  (78) 

Eq.  (78) 

Eq.  (78) 

ap/a^  =  0 

ai/a#,.  =  0 

dcojd^,  =  0 

dvjdl=0 

O 

m 

forces) 

o 

28 

Impermeable  slip  (no 

Eq.  (78) 

Eq.  (78) 

Eq.  (78) 

dP/d^,  =  0 

a/t/ai,  =  0 

dcold^,  =  0 

a^/ai,  =  0 

'c/5 

forces) 

CLh 

30 

Free  surface 

Eq.  (34) 

Eq.  (34) 

Eq.  (35) 

Eq.  (33) 

a/t/a#,  =  0 

dcojd^,  =  0 

dvjd4,=0 

40 

Zero  gradient 

etz/ei,  =  0 

er/ai,  =  o 

=  0 

dPjd^,  =  0 

ayt/ai,  =  0 

dcold^,  =  0 

dv,ld4  =  0 

41 

Translational  periodicity, 
w/  ghost  cells 

* 

* 

* 

* 

* 

* 

* 

42 

Translational  periodicity, 
w/o  ghost  cells 

* 

* 

* 

* 

* 

* 

* 

43 

Pole  (I-around) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

44 

Pole  (j-around) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

45 

Pole  (k  around) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

Eq.  (80) 

c/5 

50 

Cylindrical  zero  gradient 

* 

* 

* 

* 

* 

* 

* 

■a 

g 

51 

Rotational  periodicity,  w/ 

* 

* 

* 

* 

* 

* 

* 

m 

ghost  cells 

52 

Rotational  periodicity,  w/o 

* 

* 

* 

* 

* 

* 

* 

_o 

ghost  cells 

p 

9- 

60 

No-slip/centerplane 

* 

* 

* 

* 

* 

* 

* 

S 

o 

61 

x-axis  symmetry 

0 

dVjd^  =  0 

aiF/a^.  =  0 

ap/af,  =  0 

ai/ai,  =  0 

dajd^,  =  0 

dvjd4=0 

62 

y-axis  symmetry 

dUjd^,  =  0 

0 

dWjd^^  =  0 

dPjd^,  =  0 

ayt/ai,.  =  0 

dcold^,  =  0 

dvjd4=0 

63 

z-axis  symmetry 

dUjd^,  =  0 

er/ai,  =  o 

0 

dP/d^,  =  0 

ayt/ai,  =  0 

dco/d^,  =  0 

dvjd4=0 

91 

Multi-block  w/  ghost  cells 

* 

IN 

* 

* 

* 

* 

* 

92 

Multi-block  w/o  ghost  cells 

* 

* 

* 

* 

* 

* 

* 

99 

Blanked  out  points 

0 

0 

0 

0 

0 

0 

0 

See  text  for  detailed  description 


31 


Domain  truncation  boundaries 

For  external- flow  hydrodynamies,  an  infinite  unbounded  fluid  often  represents  the 
physieal  domain.  This  requires  that  the  eomputational  domain  be  truneated  to  a  size  that  ean  be 
eeonomieally  filled  with  grid  points,  but  has  no  influenee  on  the  eomputed  solution.  Aetual 
loeation  of  boundaries  and  influenee  on  solution  must  be  evaluated  as  part  of  the  verifieation  grid 
studies  whieh  is  diseussed  in  Section  7.  However,  BC  types  used  on  truncated  domain 
boundaries  are  listed  in  Table  7  and  represent  inlet,  exit,  far-field,  and  prescribed  boundaries. 

For  the  inlet  boundary  condition  (ibtyp=10),  the  velocity  field  is  set  by  the  input 
parameters  UINF,  VINF,  WINF,  pressure  is  zero  gradient,  and  the  turbulence  is  set  to  the  free 
stream  values  =  1.0x10  =  9.0 .  The  user  specified  freestream-velocity  unit  vector 

defined  by  UINF,  VINF,  WINF  provides  capability  to  specify  the  angle  of  attack  to  the 
vehicle. 

The  exit  boundary  condition,  ibtyp  =  11  is  derived  assuming  that  the  boundary  is  far 

d^U  d^V  d^W 

downstream  such  that  streamwise  viscous  effects  are  zero,  i.e.,  — ;r  =  — v  =  — ^  =  0.  This 

5^-  5^/ 

allows  velocity  on  boundary  to  be  calculated  using  zero_sd  from  (77)  and  all  other 
variables  extrapolated  using  zero_fd  from  (76). 

There  are  two  far-field  conditions,  ibtyp=12  and  13.  The  latter  (ibtyp=13)  specifies 
that  velocity  field  is  set  by  the  input  parameters  UINF,  VINF,  WINF  and  pressure  and 
turbulence  variables  are  zero  gradient.  This  is  preferred  option,  but  requires  that  boundary 
location  be  sufficiently  far  from  vehicle.  The  former  (ibtyp=12),  on  the  other  hand,  sets  the 
axial-component  of  velocity  to  UINF  and  pressure  to  zero  while  all  other  variables  are  assumed 
to  be  zero  gradient. 

Prescribed  boundary  condition,  ibtyp=14,  must  be  used  in  combination  with  a  user¬ 
generated  restart  file,  format  of  which  is  documented  in  Appendix  A.6.  This  condition  holds  all 
variables  constant  to  those  in  the  restart  file  except  for  pressure,  which  is  calculated  assuming  a 
zero-gradient  condition.  Typically,  this  condition  is  used  when  specifying  either  an  analytical 
(e.g.,  Blasius  fiat  plate  boundary  layer)  or  previously  computed  flow.  An  example  of  the  latter 
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may  be  an  isolated  propeller  blade  simulation,  where  the  wake  from  upstream  is  eomputed  by  a 
previous  simulation,  and  used  to  preseribe  inflow  boundary  eonditions. 

Physical  boundary  conditions 

Physieal  BC’s  are  due  to  either  solid  surfaees  or  water-air  free  surfaee.  Available  options 
are  listed  in  Table  7. 

There  are  two  no-slip  boundary  eonditions.  The  first  (ibtyp=20)  is  for  surfaees  whieh 
are  moving  with  the  grid  and  the  seeond  (ibtyp=22)  is  for  surfaees  whieh  are  moving  in  the 
relative-frame  system  provided  angular  veloeities  (agvx  or  agvz),  whieh  are  defined  in  the 
namelist  input  file  cf  d_ship  .  nml. 

Impermeable  slip  boundary  eonditions  are  speeified  using  ibtyp=27  and  28,  the  only 
differenee  being  whether  or  not  the  specified  boundary  be  included  (27)  or  not  (28)  in  the 
calculation  of  forces  and  moments.  This  distinction  is  useful  since  impermeable  slip  boundaries 
are  often  used,  for  example,  to  model  water-tunnel  walls  or  stream  surfaces  and  it  may  not  be 
desired  to  include  the  forces  on  these  boundaries  in  the  integration  process.  The  boundary 
condition  is  formulated  using  contravariant  velocity  components  at  one  point  off  the  boundary 
and  by  forcing  the  normal  component  to  be  zero.  For  example,  on  a  j=l  slip  surface,  the 
velocity  boundary  conditions  are 

U{i,\,k)  =  +  U^x^ 

V{i,\,k)  =  U^y^+U^y^+U^y^  (78) 

W{i,l,k)  =  U^z^+U^z^+U^z^ 

where 

U^=j(u  {i,  2,  k)bl  +  V{i,2,k)bl+  W{i,  2,  k)bl) 

U^=0  (79) 

U^=j(u  {i,  2,  k)b^+  V{i,  2,  k)  h'  +  W{i,  2,  k)bl) 

and  where  all  other  variables  are  calculated  using  a  zero  gradient  condition.  Similar  expressions 
can  be  dervied  for  i=constant  and  k=constant  surfaces. 

As  indicated  in  Table  7,  free  surface  (ibtyp=30)  BC  have  been  described  earlier  in 
Section  3. 
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Computational  boundaries 

Computational  boundaries  are  due  to  grid  topology,  modeling  assumptions,  and  multi- 
bloek  domain  deeomposition.  Available  options  in  CFDSHIP-IOWA  are  listed  in  Table  7  and 
inelude  zero-gradient,  translational  and  rotational  periodicity,  pole  singularities,  symmetry,  and 
multi-block  boundaries.  Each  type  is  now  described. 

Zero-gradient  boundaries  (ibtyp=4  0)  assume  that  all  variables  have  zero  gradient 
behavior.  This  condition  is  provided,  but  not  often  used  because  symmetry  conditions,  which  set 
the  normal  component  of  velocity  to  zero,  are  typically  more  appropriate. 

Translational  periodicity,  as  shown  in  Figure  8,  can  be  prescribed  using  either  ibtyp  = 
41  or  42,  the  difference  being  that  the  former  uses  2  ghost  cells  on  each  boundary  to  maintain 
solver  order-of-accuracy  and  the  latter  uses  a  simple  weighted  average  of  the  adjacent  field  point 


(a)ibtyp=42  (b)ibtyp=41 


Figure  8.  Translational  periodic  boundary  conditions,  ibtyp=41  &  42. 


values.  This  boundary  condition  can  be  used  in  either  Cartesian  or  cylindrical  coordinates.  Note 
that  all  flow  variables  are  assumed  periodic. 

Pole  boundary  conditions  are  based  upon  simple  average  of  variables  one  grid  point  off 
the  pole  in  the  interior  of  the  domain  with  three  orientations  available:  i-around  ( ibtyp=43),  j- 
around  (ibtyp=44),  and  k-around  (ibtyp=45).  The  orientation  specifies  the  summation 
index,  e.g.,  for  k-around  and  ibdir=-l-2,  the  U  component  of  velocity  on  the  pole  would  be 


34 


1 


(80) 


U{iJ,k)  = 


kbce 

-  z  U{i,j  +  \,k) 

J-  k=lrhr‘v 


kbce  -  kbcs  + 1  t=kbcs 

Rotational  periodicity  about  the  x-axis,  as  shown  in  Figure  9,  ean  be  preseribed  using 
either  ibtyp  =51  or  52.  As  with  translational  periodieity,  the  differenee  between  the  two  is  that 
the  former  uses  2  ghost  eells  on  eaeh  boundary  to  maintain  solver  order-of-aeeuraey  and  the 
latter  uses  a  simple  weighted  average  of  the  adjaeent  field  point  values.  However,  in  this  ease, 
sinee  the  radial  and  eireumferential  eomponents  are  the  periodie  variables,  V  and  W  veloeity 
eomponents  are  ealeulated  using  the  following  transformation, 

V=V  cos  d  +  W  sin  0 


Va  =-Vsm0  +  W  cos  0 


(81) 


0  =  tan  \z/y) 

Note  that  this  boundary  condition  is  restricted  to  Cartesian  grids  (icoord=l/2)  and  with 
periodicity  about  the  x-axis  only. 


Figure  9.  Rotational  periodic  boundary  conditions,  ibtyp=51  &  52. 


Centerplane  condition  (ibtyp=60)  is  a  boundary  condition  used  for  half-domain  ship- 
hull  simulations  where  part  of  the  boundary  is  no-slip  and  part  is  centerplane,  as  shown  in  Figure 
10.  Demarcation  between  the  no-slip/centerplane  is  determined  by  testing  the  y-coordinate,  i.e., 
if  |y|  <  1.0x10  then  the  boundary  point  is  assume  to  be  centerplane,  where  all  variables  are 
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Figure  10.  Centerplane/no-slip  boundary  condition,  ibtyp=60. 

calculated  using  zero_fd,  except  for  the  V-component  of  velocity  which  is  set  to  zero, 
otherwise  the  point  is  treated  as  no-slip.  This  boundary  condition  is  typically  used  to  resolve 
bow,  stem,  and/or  keel  using  a  staircase  resolution  (e.g..  Stem  et  ah,  1996).  It  is  restricted  to 
Cartesian  grids  (icoord=l/2). 

Symmetry  conditions  (ibtyp=61,  62,  63)  are  available  for  each  coordinate  direction. 
Formulation  is  simple:  normal  component  of  velocity  is  set  to  zero  (61,  U=0;  62,  V=0,  63, 
W=0),  and  all  other  flow  variables  are  assumed  to  be  zero  gradient. 

CFDSHIP-IOWA  utilizes  2  types  of  multi-block  boundary  conditions,  pointwise 
continuous  (i.e.,  abutted-block  interface)  and  overset  (i.e..  Chimera),  examples  of  both  are  shown 
in  Figure  11. 
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Figure  1 1 .  Overset  and  patehed  multiblock  grids  for  airfoil. 


For  abutted  multi-block  boundary  conditions,  two  options  are  available,  ibtyp=91  and 
92,  where  the  former  includes  2  rows  of  ghost  cells,  which  are  added  after  the  grid  is  read  from 
file.  For  ibtyp=91,  all  independent  flow  variables  in  the  ghost  cells  are  computed  as  field 
variables  in  the  donor  block.  This  conserves  mass  and  momentum  across  the  boundary. 
Furthermore,  by  using  2  rows  of  ghost  cells,  the  5-point  stencil  and  therefore  the  solver  order-of- 
accuracy  is  maintained  on  the  boundary.  In  contrast,  ibtyp  =92  calculates  boundary  values 
using  a  weighted  average,  based  upon  distance  from  boundary,  of  the  field  values  on  each  side  of 
the  interface.  This  condition  is  essentially  a  first  order  treatment  across  the  interface  and  does 
not  solve  governing  equations  for  these  boundaries.  It  is  noted  that  a  consistent  approach  must 
be  used,  i.e.,  one  or  the  other  must  be  used  for  all  abutted  boundaries,  mixed  usage  is  not 
currently  supported.  Moreover,  for  complex  grid  topologies,  it  is  possible  that  ibtyp  =91  may 
fail  to  properly  identify  ghost  cells.  In  those  cases,  ibtyp  =92  must  be  used. 

For  overset  multi-block  boundary  conditions.  PEGASUS  software  from  NASA  Ames 
Research  Center  must  be  used.  Interface  and  implementation  are  described  in  the  following  sub¬ 
section. 
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4. 7  Chimera  overset  gridding 

Capability  for  simulations  using  Chimera-style  overset  domain  deeomposition  is 
aeeomplished  through  interfaee  with  PEGASUS  version  5.1  (Suhs,  et  ah,  2001),  whieh  is  the 
latest  version  of  the  PEGASUS  series  of  mesh  interpolation  eodes,  originally  developed  at 
NASA  Ames  Researeh  Center.  The  main  purpose  for  the  development  of  version  5.1  was  to 
deerease  the  number  of  user  inputs  required  and  to  allow  for  easier  operation  of  the  eode.  A 
basie  deseription  of  Chimera  methodology  is  deseribed  in  the  Version  4  manual  (Suhs  and 
Tramel,  1991).  It  should  be  noted  that  PEGASUS  is  restrieted  to  U.S.  institutions  and 
researehers  due  to  export  eontrol  regulations.  Other  options  for  eomputation  of  interpolation 
eoeffieients  inelude  CHALMESH  (Petersson,  1999),  OVERTURE  (Quinlan  et  ah,  2002)  and 
PEGISUE  (Denny,  2002). 

CEDSHIP-IOWA  is  designed  to  use  double-fringe  hole  and  outer  boundaries  and  level-2 
interpolation.  Eigure  12  shows  an  example  of  an  overset  grid  system  for  an  airfoil  to  aid  in 
defining  terminology.  Comparing  Eigures  12a  and  12b,  it  ean  be  seen  that  double-fringe  blanks 
out  2  layers  of  eells  around  the  outer  and  hole  boundaries.  Both  layers  are  interpolated  from  the 
donor  mesh.  This  permits  use  of  the  normal  5-point  steneil  for  all  field  points  and  maintains 
order  of  aeeuraey  near  boundaries.  Single-fringe  boundaries,  on  the  other  hand,  require  use  of  3- 
point  steneils  for  the  field  points  adjaeent  to  hole  and  outer  boundaries,  whieh  results  in 
reduetion  of  the  order-of-aeeuraey.  The  downside  in  using  a  double  fringe  is  that  it  requires 
more  mesh  points  and  makes  it  more  diffieult  to  obtain  the  required  overlap  between  hole 
boundaries  and  outer  boundaries. 

Eigure  12e  shows  the  impaet  of  enabling  level-2  interpolation.  This  is  a  method  where 
onee  the  minimum  hole  has  been  established;  the  hole  is  enlarged  to  improve  the  eommunieation 
among  meshes.  In  addition,  it  ean  be  seen  that  the  outer  boundary  of  the  foil  mesh  has  been 
moved  inward.  This  is  aeeomplished  in  PEGASUS  by  eomparing  relative  grid  quality  between 
interpolation  and  donor  meshes  so  as  to  reduee  mesh  disparity,  whieh  ean  severely  impaet 
aeeuraey.  In  addition,  level-2  interpolation  provides  eapability  to  ereate  holes  with  refinement 
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(a)  single  fringe  outer  and  hole  boundaries 


rsi 


□.i 

02 

a 

-□2 

-□.i 


■016 

,  X 

(c)  double  fringe  outer  and  hole  boundaries 
with  leMel-2  interpolation 


(b)  double  fringe  outer  and  hole  boundaries 


Figure  12.  Illustration  of  PEGASUS  terminology. 


meshes  that  are  added  to  the  domain.  This  eannot  be  accomplished  through  the  usual  hole¬ 
cutting  methods,  since  refinement  meshes  typically  do  not  have  solid  walls  that  would  be  used 
for  definition  of  hole-cutting  boundaries.  Refinement  meshes  have  been  used  in  the  aerospace 
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community  (e.g.,  Meakin,  R.,  1999)  for  flow  adaptive  refinement  within  eontext  of  struetured 
flow  solvers. 

Given  a  valid  PEGASUS  interpolation  file  (i.e.,  in  PEGASUS  standard  naming 
convention,  X INTOUT),  implementation  in  CEDSHIP-IOWA  is  straightforward.  This  file 
should  be  renamed  FNAMEI  .x intout  sinee  CEDSHlP-lOWA  automatieally  looks  for  this  file 
in  the  execution  directory.  If  it  exists,  the  file  is  read  and  a  summary  is  printed  in  the  standard 
output.  Otherwise,  the  code  continues  and  prints  a  notification  message  that 
FNAMEI .  xintout  does  not  exist.  It  should  be  emphasized  that  in  the  boundary  condition 
input  file  FNAMEI .  bcs  Chimera  outer  boundaries  remain  unspecified,  which  for  most 
problems,  greatly  simplifies  and  shortens  this  input  file. 

As  a  summary,  Eigure  13  provides  a  flowehart  summarizing  the  overset-grid  proeess. 
Eurther  details  can  be  found  in  the  appropriate  domain  connectivity  software,  which  in  this  case 
is  PEGASUS  version  5.1. 


4.8  Calculation  of  forces  and  moments 


The  fluid  forees  and  moments  aeting  on  all  solid  surfaees  are  obtained  by  integration  of 
the  normal  and  tangential  stresses  over  all  no-slip  (ibtyp=20,22)  and  some  of  the  user- 
identified  slip  boundaries  (i.e.,  ibtyp=27).  The  fluid  stress  tensor  on  a  no-slip  surface  is 
eomprised  of  eomponents  due  to  pressure  and  viscous  stress 


dU,  dU. 


^dXj 


-  +  - 


dx, 


(82) 


On  the  solid  surfaees,  the  fluid  forees  and  moments  are  determined  through  integration  of  (82) 

F  =  (83) 


M  =  rxF 


(84) 


where  rf  is  the  unit-vector  normal  to  a  ^-coordinate  surfaee,  dS^  is  the  loeal  surface  area 
element,  and  r  is  the  position  vector 
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(85) 


b^i  +  blj  +  bik 


l(b;)  +(«)  +(«) 
ds‘  =  ^{b;f +{bif +{bif d(’d4‘ 


and  where  are  eye  lie.  Integration  of  wetted  surfaee  area  S  and  equations  (83)-(84)  is 

aoeomplished  using  a  eell-eentered  2D  trapezoidal  rule.  Caleulation  of  the  foree  is  broken  down 
into  eontributions  from  piezometrie  pressure  (eppiezo),  hydrostatie  pressure  (ephydro),  whieh  is 
valid  only  for  Fr^tQ,  and  viscous  shear  stress  (cf).  These  contributions,  along  with  the  total  force 
(ctot),  are  written  for  each  base  coordinate  direction  every  iteration.  The  moments  are  treated  in 
a  similar  fashion. 

For  overset  grid-systems  with  overlapping  surface  meshes  on  solid  surfaces,  forces  and 
moments  computed  by  CFDSHIP-IOWA  will  contain  errors  due  to  multiply  defined  values  and 
hole  boundaries.  Special  tools  must  be  used  to  make  flow  variables  single -valued.  Currently, 
the  force  and  moment  computation  (FOMOCO)  tools  from  NASA  (Chan  and  Buning,  1996) 
represent  the  only  option  for  performing  this  task.  Two  programs  make  up  this  toolset,  MIXSUR 
and  OVERINT  where  the  former  creates  a  single-valued  mixed  surface  with  triangular  zipper 
grids  and  the  latter  performs  integration  over  this  new  surface.  Unlike  the  NASA  flow  solver 
OVERFLOW,  which  reads  MIXSUR  output  as  in  input  file  and  computes  forces  and  moments 
on  the  fly  by  calling  OVERINT  as  a  subroutine,  CEDSHIP-IOWA  currently  uses  these  tools  in  a 
stand-alone  post-processing  fashion.  Euture  versions  may  more  tightly  integrate  with  these  tools 
or  appropriate  alternatives. 
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4.9  Algebraic  equation  solver 

The  overall  method  is  fully  implicit  and  there  are  four  locations  in  the  code  that  require 
iterative  solvers:  momentum  predictor  step  (59);  pressure  equation  (64);  turbulence  model 
equations  (24);  and  KFSBC  (31).  Currently,  a  line- ADI  scheme  with  a  pentadiagonal  solver  and 
under  relaxation  is  used  to  solve  the  algebraic  equations.  The  pentadiagonal  solver  is  modified 
to  account  for  hole  boundaries 

а,  =  apiBLANK, 

б.  =b,niBLANK. 
c.=c,Jl\-IBLANK,) 

^  ’  (86) 

=  d^niBLANK. 

e.  =  e.nBLANK, 

f  =  f^BLANK^+{\-lBLANK,)AI), 

where  IB  LANK  is  0  for  a  hole  or  fringe  point  and  1  for  a  field  point,  ai,  bi,  Cj,  dj,  e,  correspond  to 
5  diagonals  of  the  pentadiagnoal  matrix,  ft  is  the  right  hand  side,  and  corresponds  to  the  flow 
variable  on  the  hole  or  outer  boundary  fringe. 

4.10  Algorithm  summary  and  flowchart 

Figures  14  and  15  summarize  the  algorithmic  structure  of  CFDSHIP-IOWA  in 
flowcharts,  including  pre-  and  post-processing  functions,  major  time  marching  and  iteration 
loops,  and  input  and  output  files. 
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Figure  15.  Detailed  flowchart  for  PISO  velocity-pressure  coupling  solver 
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5.  CODE  DEVELOPMENT  AND  HIGH-PERFORMANCE  COMPUTING 

Development  of  CFDSHIP-IOWA  has  taken  plaee  over  the  past  8  years,  a  period  in 
whieh  high-performanee  eomputing  (HPC)  platforms  have  evolved  from  Cray  YMP  vector 
processors  to  the  array  of  architectures  found  today  which  includes  commercial  off-the-shelf 
(COTS)  Beowulf  clusters,  SGI  scalable  distributed-shared  non-uniform  memory  access  (NUMA) 
architectures,  and  IBM  pipelined  superscalar  architectures.  Developing  application  codes  which 
are  portable  and  which  are  capable  of  harnessing  a  given  platform’s  capability,  is  a  challenging 
task.  In  addition,  CFDSHIP-IOWA  had  to  meet  other  objectives  such  as  supporting  student 
theses  and  project  research,  as  well  as  vertical  transition  to  other  universities,  industry,  and 
government  labs. 

The  approach  used  to  meet  these  objectives  has  been  based  upon  a  flexible  data  structure, 
capability  for  both  serial  and  parallel  computing,  adherence  to  standards  such  as  MPI,  modem 
programming  using  Fortran90/95,  and  interface  with  existing  3'^‘^-party  software  for  grid 
generation  and  boundary  condition  specification  (GRIDGEN),  Chimera  overset  grid  methods 
(PEGASUS,  Chimera  Grid  Tools,  EOMOCO),  and  post-processing  and  visualization 
(TECPEOT).  It  is  noted  that  development  has  taken  place  concurrently  with  several  DOD 
HPCMO  Challenge  Projects  from  1996  to  the  present  (Rood,  1997;  Rood,  1998;  Rood,  1999; 
Rood,  2000;  Kim,  2001).  The  following  sections  present  an  overview  of  the  code  and  data 
stmctures,  parallel  computing,  portability,  and  distribution,  extraction,  compilation  and 
execution. 

5.1.  Code  and  data  structures 

Code  and  data  stmctures  are  critical  to  successful  implementation  of  scalable  parallel 
computing  and  portability.  Unfortunately,  highly  efficient  code  is  often  hard  to  understand  by 
users.  Because  CEDSHIP-IOWA  was  designed  to  support  research  and  development  of  new 
models,  it  is  well  documented  with  in-code  comment  statements  and  the  stmcture  at  the 
subroutine  level  is  based  upon  3D  index-ordered  arrays,  which  are  intuitively  easier  to 
understand  in  context  of  stmctured-grid  CED,  and  which  lends  itself  naturally  to  distributed 
multi-block  computing. 

Top-level  stmcture  of  the  code,  which  is  also  shown  in  Eigure  14,  is  as  follows 
•  program  cfdship  iowa 
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o  pre_process 
o  free_surface 
o  grid_transformation 
o  eddyvs 
o  piso 

o  post_processl 
o  post_process2 

In  process  of  writing  code,  several  guidelines  were  adopted.  At  top-level,  all  arrays  are  dynamic, 
ID,  and  allocatable  to  fit  current  grid  size.  At  subroutine  level,  computation  is  performed  on  a 
single  block;  all  arrays  are  3D  and  either  explicit-shape  dummy  arrays  (both  array  and 
dimensions  are  in  argument  list)  or  automatic  arrays  (arrays  which  are  not  in  argument  list  and 
which  are  created/destroyed  upon  entry/exit  of  a  given  procedure).  Index  pointers  are  used  to 
indicate  relative  location  in  ID  array  and  are  defined  as  follows. 


first ( 1 ) =1 

length  (1)  =imax  (1)  *  jmax  (1)  *k;max  (1) 
last ( 1 ) =length ( 1 ) 
ntot=length ( 1 ) 
do  m=2,nmesh 

(87) 

first (m)=last (m-1 ) +1 
length (m) =imax (m) * jmax (m) *kmax (m) 
last (m) =f irst (m) +length (m) -1 
ntot=ntot+length (m) 

enddo 

Fortran90  module  procedures  are  used  in  preference  over  COMMON  blocks  to  share  definitions 
and  values  of  data  between  program  units  due  to  ease  of  code  maintenance  and  data-hiding 
capability.  To  reduce  coding  errors,  IMPLICIT  NONE  statements  are  used  in  all  routines. 

By  design,  this  approach  easily  permits  parallel  multi-block  implementation  via  message¬ 
passing  interface  (MPI)  and  the  single-program  multiple  data  (SPMD)  paradigm.  As  shown  in 
the  code  fragment  in  Figure  16,  serial  implementation  for  a  given  function  performs  a  DO  loop 
over  each  block  with  the  locn  pointer  set  to  first  (m) .  This  passes  each  block  to  the 
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c 

c  --  calculate  transformation  metrics 

c 

#ifdef  SERIAL 

do  m  =  l,nmesh 

locn  =  first (m) 

#endif 

#ifdef  PARALLEL 

m  =  mymesh 
locn  =  1 

#endif 

call  metric ( imax (m) ,  jmax (m) ,  kmax (m) , 
xp (locn) , yp (locn) , zp (locn) , 
xpO (locn) , ypO (locn) , zpO (locn) , 
xpOO (locn) , ypOO (locn) , zpOO (locn) , 
bll (locn) ,bl2 (locn) ,bl3 (locn) , 
b21 (locn) ,b22 (locn) ,b23 (locn) , 
b31 (locn) ,b32 (locn) ,b33 (locn) , 
all (locn) , a22 (locn) , a33 (locn) , 
al2 (locn) ,al3(locn) ,a23(locn) ,aji (locn) , 
fl (locn)  , f2 (locn) , f3 (locn) , 
agvx (m) , agvy (m) , agvz (m) , 
dxdt (locn) , dydt (locn) , dzdt (locn) ) 

#ifdef  SERIAL 
enddo 

#endif 

Figure  16.  Code  fragment  illustrating  use  of  pointer  and  CPP  statements. 


subroutine  metric (imax, jmax, kmax, x,y,z,x0,y0,z0,x00,y00,z00, 
bll,bl2,bl3,b21,b22,b23,b31,b32,b33, 
all,a22,a33,al2,al3,a23,aji,fl,f2,f3, 
agvx,  agvy, agvz, dxdt, dydt, dzdt) 
use  global_parameters 
implicit  NONE 

integer  i, j , k, imax, jmax, kmax 

real  (kind=double) ,  dimension (imax, jmax, kimax) ; :  x, y , z, xO, yO, zO , xOO , yOO , zOO 

Figure  17.  Subroutine  fragment  illustrating  3D  explicit-shape  dummy  arrays. 

subroutine,  as  shown  in  Figure  17,  in  a  sequential  fashion.  It  is  noted  that  grid  and  flow  variable 
arrays  are  dimensioned  at  the  top-level  to  hold  the  entire  system,  i.e.,  to  ntot .  In  contrast,  for 
parallel  implementation,  the  DO-LOOP  is  eliminated  and  the  pointer  FIRST  is  set  to  1  since 
each  processor  only  has  the  data  of  a  single  block.  As  such,  each  processor  executes  its  own 
copy  of  the  executable.  Except  for  input  and  output  routines,  communication  between 
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processors  occurs  only  at  the  boundary  condition  subroutines.  Processor  0,  through  the  use  of 
branching  statements,  handles  control  of  input  and  output.  It  should  be  noted  that  the  suffix  (.F 
or  .F90)  on  the  source  code  files  indicates  that  the  file  contains  CPP  statements  which  in  turn 
indicates  the  file  contains  parallel-specific  programming  statements. 

5.2  Parallel  Computing 

CFDSHIP-IOWA  achieves  scalable  parallel  performance  using  several  parallel  models. 
Originally,  it  was  designed  as  a  distributed-memory  coarse-grain  message-passing  model  based 
upon  domain  decomposition  and  the  message-passing  interface  (MPI).  For  good  performance, 
this  approach  requires  static  load  balancing,  i.e.,  the  grid  system  be  decomposed  into  nearly 
equal  sized  blocks,  and,  in  addition,  requires  a  block  for  each  processor.  While  very  efficient, 
this  process  becomes  tedious  and  makes  post-processing  difficult,  especially  when  large  numbers 
of  processors  (>32)  are  required.  To  alleviate  this  problem,  a  second  parallel  model  using 
OpenMP  threads  for  shared-memory  fine-grain  (i.e.,  loop-level)  parallelism  was  introduced. 
This  model  is  used  in  combination  with  MPI  to  achieve  multi-level  parallelism,  which  permits 
use  of  very  large  numbers  of  processors  and  can  perform  dynamic  load  balancing.  In  the 
following,  each  of  these  models  is  discussed  in  detail. 

To  demonstrate  the  performance  of  the  message-passing  model,  a  105x61x30  single¬ 
block  grid  for  the  Wigley  hull  is  decomposed  into  2-,  4-,  8-,  12-,  16-,  and  32-block  grid  systems 
and  simulations  timed  using  both  ibtyp=91  and  ibtyp=92  (i.e.,  with  and  without  ghost  cells, 
respectively)  multi-block  boundary  conditions.  Figure  18  shows  the  impact  of  decomposition  on 
both  the  total  problem  and  individual  block  sizes.  Both  boundary  conditions  impose  an  overhead 
due  to  replicated  points,  but  overhead  when  using  ibtyp=91  increases  significantly  with 
processors  due  to  increasingly  large  number  of  ghost  cells.  This  directly  impacts  the  scalability 
of  memory  utilization.  However,  this  tradeoff  is  often  necessary  if  the  improved  accuracy  of 
ibtyp=91  is  required. 

Figure  19  shows  the  parallel  speedup  =1^/7),  (n),  where  Ts  is  the  single  processor 

wall-clock  time  and  Tp(n)  is  the  wall-clock  time  for  n  processors,  for  both  multi-block  boundary 
conditions  on  both  the  CRAY  T3E  and  the  SGI  Origin  2000  (02K).  For  ibtyp=9I, 
degradation  of  speedup,  due  to  above-mentioned  overhead,  is  obvious,  especially  for  16  and  32 
processors.  Agreement  between  machines  for  ibtyp=92  is  not  consistent,  i.e.,  the  02K  shows 
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almost  linear  speedup  through  32  processors  and  the  T3E  shows  some  drop  off  for  n  =  16  and  32. 
Since  timings  were  not  made  in  dedicated  mode,  the  reproducibility  is  dependent  upon  system 
load,  which  may  explain  the  T3E  performance.  Einally,  it  is  noted  that  the  timings  show  that 
99.5%  of  the  code  is  parallelized. 


Eigure  18.  Impact  of  domain  decomposition  and  use  of  ghost  cells  on  number  of  grid  points. 


(a)  CRAYT3E 

Eigure  19.  Parallel  speed-up:  WigleyHull. 
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(b)  SGI  Origin  2000 


Shared-memory  fine-grain  (i.e.,  loop  level)  parallelism  is  introdueed  using  OpenMP 
whieh  was  designed  to  exploit  eertain  eharaeteristies  of  shared-memory  arehitectures.  Systems 
that  don't  fit  the  classie  shared-memory  architeeture  (e.g.,  Beowulf  elusters)  may  use  OpenMP 
but  typical  performance  is  very  poor  due  to  high  latencies  in  communication.  Currently, 
CFDSHIP-IOWA  implementation  of  OpenMP  is  only  supported  on  the  SGI  Origin  and  IBM  SP 
machines  through  the  use  of  the  automatic  parallelization  option  (-apo)  in  the  Fortran  90 
compiler.  In  particular,  the  SGI  Origin  uses  a  distributed-shared  memory  (DSM)  architecture, 
which  can  effectively  utilize  the  shared-memory  model. 

To  demonstrate  performance,  parallel  simulations  were  performed  on  the  NAVO  02K 
using  a  surface-piercing  fiat  plate  (SPFP)  with  grid  dimensions  of  105x61x30  and  with  number 
of  processors  ranging  from  1  to  32.  Figure  20  shows  that  speedup  stalls  at  about  4  for  this 
problem  on  this  machine.  Analysis  of  cache  utilization  shows  that  the  code  is  “non-cache 
friendly”  for  scalable  shared-memory  applications.  While  this  is  an  area  of  work  for  the  future, 
even  a  speed-up  of  4  is  useful  for  dynamic  load  balancing. 


Figure  20.  Parallel  speedup,  shared  memory  on  NAVO  02K. 


Mixed-mode,  or  multi-level,  parallelism  utilizes  both  MPI  and  OpenMP  to  achieve 
dynamic  load  balancing.  Relative  block  size  and  the  total  number  of  processors  specified  by  the 


51 


input  variable  total_num_procs  sets  the  number  of  OpenMP  threads  that  are  used  for  each 
block, 


num  thrd  =  max  <  int 


float  {ntot) 


\ 


total  _  num  _  procs 
^  float  [ntot  _  sum)  j 


,1 


(88) 


Table  8  shows  an  example  for  a  4-block  grid  system  for  an  open-water  propeller.  The  largest 
block  is  8.5  times  bigger  than  the  smallest  block  and  comprises  65%  of  the  total  number  of  grid 
points.  The  distribution  of  processors  for  total_num_procs  =  16  and  32  are  shown. 
Unfortunately,  performance  studies  have  not  yet  been  undertaken  for  a  simple  single-block  grid 
such  as  the  SPFP  or  the  Wigley  Hull.  However,  based  upon  the  OpenMP  results,  it  should  be 
expected  that  using  10  or  20  processors  on  Block  4  may  be  inefficient  and  that  domain 
decomposition  should  be  used  to  bring  it’s  sub-blocks  closer  to  the  size  of  blocks  1-3. 


Table  8.  Dynamic  load  balancing  for  4-block  grid  system. 


block 

block  size  (ntot) 

numthrd 

total_num_procs  =  16 

numthrd 

total_num_procs  =  32 

1 

102951 

2 

4 

2 

102951 

2 

4 

3 

56203 

1 

2 

4 

480751 

10 

20 

ntotsum  =  742856 

Total  procs  used  =  15 

Total  procs  used  =  30 

5.3  Portability 

The  code  was  designed  to  be  portable  across  the  range  of  machines  currently  available  at 
the  DOD  High-Performance  Computing  Modernization  Program  (HPCMP)  centers.  Currently 
this  includes  the  SGI  Origin  2000  &  Origin  3000,  Cray  T3E,  Cray  SVl,  Cray  T90,  and  IBM  SP2. 
In  addition,  the  code  has  been  compiled  on  the  DEC  Alpha,  HP  Workstation,  and  Intel  X86 
personal  computers,  the  latter  of  which  uses  Compaq  Visual  Eortran. 

Portability  is  achieved  through  the  use  of  MPI,  the  C  preprocessor  CPP,  and  the  UNIX 
make  utility.  Vendor-optimized  MPI  is  available  on  nearly  all  platforms.  Otherwise,  the 
generic  MPICH  libraries  can  be  used.  CPP  is  used  to  build  either  the  parallel  or  serial  versions 
from  a  single  source  code.  MPI-  and  serial-specific  code  is  isolated  through  the  use  of 
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'#IFDEF  SERIAL’,  '#IFDEF  PARALLEL’,  and  ' #END IF '  CP P  directives.  Although  not 
typically  required,  the  source  code  for  either  version  can  be  readily  obtained  by  simply  running 
CPP  with  the  appropriate  directives  (i.e.,  -DSERIAL  or  -DPARALLEL).  In  addition,  directives 
are  used  to  set  MPI  data  types  to  MP  I_REAL  for  all  machines  except  for  the  CRAY  T3E,  which 
requires  MPI_DOUBLE_PRECISION.  The  code  is  written  in  FORTRAN  and  compiles  with 
Fortran  90.  The  makefile  will  build  platform  specific  versions  by  invoking  the  correct  compiler 
options  and  the  CPP  directives.  The  machines  and  options  available  are  shown  in  Table  9. 


Table  9.  Supported  platforms  and  make  options. 


Machine 

make  argument  (Serial/Parallel) 

SGI  Origin  2000 

02K/02K  MPI 

SGI  Power  Challenge  Array 

PCA/PCA  MPI 

CRAY  T3F 

T3F/T3E  MPI 

CRAY  T90 

T90/NA 

DFC  Alpha 

DEC/DEC  MPI 

HP  workstation 

HP/NA 

Intel  X86  Processor 

None  (Visual  Fortran) 

Beowulf  Cluster  (Portland  Group  Compiler) 

CEUSTER/CEUSTER  MPI 

5.4  Code  distribution,  extraction,  compilation  and  execution 

CFDSHIP-IOWA  source  code  is  distributed  as  a  compressed  tar  fide,  which  can  be 
uncompressed  and  extracted  using  the  following  UNIX  command, 

%tar  -xvf  I  uncompress  cfdship . tar . Z 

The  following  files  are  included  in  the  distribution 


filename 

Description 

readme 

Description  of  files 

makefile 

UNIX  makefile 

cfdship  mods.F90 

Definition  of  Fortran  modules 

cfdship  chimera. F90 

Interface  with  PEGASUS 

cfdship  bcs.f 

Boundary  conditions 

cfdship  frsf.f 

Free  surface 

cfdship  ke . F 

Discretization/solution  of  k-s  turbulence  model 
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cfdship  ko . F 

cf  dship_inain .  F 

cfdship  mom.f 

cfdship  mpsr.F 

cfdship  pres.f 

cf dship_ship . f 

cfdship  stio.F 

cfdship  turb.F 

cfdship  util.f 

Tools/ convert_chimera . F90 

Tools /decomp . f 

Tools/grd-atob . f 90 

Tools/grid-btoa . f 90 

Tools/ grid- re strict.! 

Tools/inlet-rst.f 

Too Is /pro long. f 

Tools/redist .  f 90 
Tools/slicer.f90 
Tools/ xyz-to-xrt .  f 


Discretization/solution  of  k-o)  turbulence  model 
Main  program 

Discretization/solution  of  momentum  equation 
Multi-processor  versions  of  subroutines 
Discretization/solution  of  pressure  equation 
Propeller  body  force,  forces  and  moments 
Setup  and  input/output  routines 
Turbulence  closure 
Solvers,  interpolation  routines,utilities 

Converts  XINTOUT  file  to  either  ASCII  or  unformatted  format 
Partitions  grid  into  sub-blocks  for  parallel  computing 
Converts  grid  file  from  ASCII  to  unformatted  format 
Converts  grid  file  from  unformatted  to  ASCII  format 
Creates  coarse  and  medium  grids  given  fine  grid 
Creates  restart  file  given  inlet  profile 

Interpolates  coarse-grid  solution  onto  medium  and  fine  grids  and 

creates  restart  files 

Changes  distribution  of  grid  points 

Extracts  subset  of  data  from  restart  file 

Convert  grid  from  Cartesian  to  cylindrical  coordinates 


The  file  extension  indicates  the  presence  of  CPP  statements  (i.e.,  .F  or  .F90)  and  whether  the 
code  is  written  in  fixed  (i.e.,  .f  or  .F)  or  free-formats  (i.e.,  .f90  or  .F90). 

After  uncompressing  and  extracting  all  files,  the  flow  code  is  compiled  using  the 
makefile  along  with  machine  specific  option  listed  in  Table  9.  For  example,  to  compile  the  flow 
code  on  the  SGI  Origin  2000  with  MPI  parallelism,  the  following  should  be  typed 

%make  02K  MPI 

Execution  is  platform  dependent  and  requires  site-specific  preparation  of  job  scripts  and 
queue  submission.  However,  serial  jobs  are  executed  with  the  following  command, 

%cfd_ship  >standard_output_file 

whereas  parallel  jobs  are  executed  using  the  mpirun  command 
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%mpirun  -np  12  cfd  ship  >standard_output_file 

In  either  mode,  the  code  automatically  looks  for  the  NAMELIST  input  file  cfd_ship .  nml, 
preparation  and  contents  of  which  are  described  in  the  next  section. 

6.  CREATING  INPUT  FILES  AND  POST-PROCESSING 

The  fourth  and  sixth  steps  of  the  CFD  Process,  as  defined  in  Section  2,  are  focused  on 
creating  input  files  and  post-processing  of  simulation  results,  respectively.  While  these  steps  are 
common  to  all  CFD  codes,  the  mechanics  of  undertaking  these  tasks  are  usually  unique  and  can 
potentially  represent  a  significant  amount  of  time  in  the  overall  CFD  Process.  Because  of  this, 
an  overall  framework  and  graphical  user  interface,  most  notably  in  commercial  codes,  is  often 
utilized  to  simplify  the  process  and  hide  low-level  details  from  user.  Another  approach  is  to  use 
scripting  languages  to  automate  repetitive  tasks  like  grid  generation  and  setting  of  parameters.  In 
recognition  of  the  challenge  this  presents  to  small  groups  developing  research  codes,  there  are 
current  efforts  sponsored  by  HPCMP  CHS  SI  to  develop  extensible  frameworks  for  managing 
input  file  creation,  file  formats,  and  communication  between  tools.  Unfortunately,  until  such 
software  is  available  to  the  CFD  community,  users  of  CFDSHIP-IOWA  will  be  required  to  have 
first  hand  knowledge  of  how  to  create  and  manipulate  data  files. 

6.1  Input  files 

CFDSHIP-IOWA  reads  5  different  input  files  to  provide  data  for  computational  grid, 
initial  conditions,  boundary  conditions,  flow  conditions,  selection  of  models,  and  specification  of 
numerical  parameters  and  post-processing  variables.  Table  10  summarizes  the  filenames, 
Fortran  unit  numbers,  and  descriptions  of  the  input  files.  Detailed  presentation  of  the  file  format 
is  provided  in  Appendix  A.  Here,  the  purpose  of,  and  the  method  for  creating,  each  file  is 
described. 
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Table  10.  Input  files 


Filename 

Unit  Number 

Description 

cfd_ship .  nml  (fixed  filename) 

8 

Namelist  input  for  runtime  variables 

grid  file  (no  filename  restrictions) 

15 

Grid  file  in  ASCII  Plot3D  format 

FNAMEO.bcs 

25 

Boundary  condition  data 

FNAMEO . xintout 

45 

Chimera  interpolation  coefficients 

FNAMEI.rsto 

35 

Restart  file  from  previous  simulation 

The  master  control  file  for  CFDSHIP-IOWA  is  named  cf  d_ship  .  nml.  Its  purpose  is 
to  specify  initial  conditions,  flow  conditions,  selection  of  models,  and  specification  of  numerical 
parameters  and  post-processing  variables.  All  input  parameters  belong  to  one  of  nine  namelists 
(control,  f low_parameters ,  grid  parameters,  iteration,  solver, 
turbulence,  f ree_surf ace,  propeller,  filenames).  A  namelist  is  a  list  of 
variable  names  that  are  always  read  or  written  as  a  group.  Namelist  input  provides  a  convenient 
interface  for  a  research  CFD  code  since  default  values  can  be  set  for  most  variables,  and  new 
variables  can  be  added  in  future  versions  without  making  previous  input  files  obsolete.  In 
general,  cf  d_ship  .  nml  is  created  by  copying  an  existing  file  and  modifying  variables  using  a 
text  editor  as  required.  Detailed  description  of  the  namelists,  variables,  and  default  values  are 
included  in  Appendix  A.2  and  an  example  can  be  found  in  Appendix  B.l. 

The  grid  file  contains  (x,  y,  z)  or  (x,  r,  6)  coordinates  of  the  structured-grid  multi-block 
system.  Note  that  for  Cartesian  and  cylindrical-polar  grids,  the  icoord  variable  in  namelist 
CONTROL  must  be  set  to  (1/2)  or  (3/4),  respectively.  The  grid  file  format  is  ASCII  PlotSD,  the 
details  of  which  are  specified  in  Appendix  A.l.  Method  for  generating  grid  file  is  at  the 
discretion  of  the  user. 

The  boundary  condition  file  specifies  boundary  condition  types  on  all  faces,  which  may 
be  arbitrarily  broken  into  rectangular  sub-patches,  of  the  computational  domain,  including  multi¬ 
block  interfaces.  For  each  patch,  the  following  information  must  be  specified  in  the 
FNAMEI .  bcs  input  file:  ibtyp  is  the  boundary  condition  type,  ibdir  is  the  inward  pointing 
normal  direction  in  computation  coordinate  direction  (+1/-1,  +11-1,  or  +3/-3 )  ,  ibcs ,  ibce , 
jbcs,  jbce,  kbcs,  kb ce  are  the  starting  and  ending  indices  in  the  coordinates 

directions,  ibcord  is  a  flag  used  to  set  the  discretization  order-of-accuracy  for  Neumann 
boundary  conditions,  i.e.,  ibcord=0  for  first  order  and  1  for  second  order,  and  if  sf  liter  is 


56 


a  flag  which  is  only  used  for  free-surface  boundaries  and  sets  the  filter  type  as  deseribed  in 
Seetion  4.4.  For  multi-block  and  periodic  conditions,  the  following  additional  information  is 
required:  ndmesh  is  the  donor  bloek  number,  idbdir  is  the  inward-pointing  normal  in  the 
donor  bloek,  and  ides,  idee,  jdes,  jdee,  kdes,  kdee  are  the  starting  and  ending 
indiees  for  the  donor  bloek.  Detailed  deseription  of  the  file  format  is  provided  in  Appendix  A.3 
and  an  example  is  provided  in  Appendix  B.3.  The  reeommended  proeedure  for  setting  boundary 
eonditions  and  ereating  FNAMEI .  bes  is  to  use  GRIDGEN  software  from  Pointwise,  Ine. 
wherein  CFDSHIP-IOWA  is  one  of  the  supported  analysis  software  options  (AS/W).  This 
allows  use  of  the  GRIDGEN  graphieal  user  interface  which  greatly  reduces  time  and  errors. 

The  final  two  input  files  are  optional.  Restart  files  are  only  required  if  initial  eonditions, 
or  preseribed  boundary  eonditions  (ibt  yp=l  4),  are  set  by  previous  simulations.  If  a  restart  file 
is  to  be  read,  the  variable  mode  in  NAMEEIST  CONTROL  must  be  set  to  1.  Restart  file  is 
typically  created  by  previous  simulation,  however,  users  can  write  eustom  software  similar  to  the 
provided  tool  inlet-rst.f  which  can  write  a  restart  file  for  speeifying  an  inlet  profile.  The 
file  eontaining  Chimera  overset  grid  interpolation  eoeffieients  is  required  only  when  using 
overset  grids.  This  fide  is  created  following  flowchart  depicted  in  Figure  13. 

6.2  Output  files  &  post-processing 

Output  fdes  provide  aeeess  to  the  simulation  results  and  ean  be  used  to  assess  iterative 
eonvergenee,  determine  forees  and  moments,  and  analyze  details  of  the  flow  field.  As  shown  in 
Table  11,  there  are  8  output  files  that  eontain  simulation  results,  however,  the  restart  file  is 
typieally  not  used  for  post-proeessing.  Four  of  the  files  are  used  for  assessing  iterative 
eonvergenee.  FNAMEO .  res  eontains  residuals,  average  divergenee,  and  evaluation  of  mass  and 
momentum  balanee  over  the  eomputational  domain,  using  the  Reynolds-transport  theorem,  at 
eaeh  time  step  (or  global  iteration).  FNAMEO .  forces  and  FNAMEO .  moments  provide  forees 
and  moments  aeting  on  the  no-slip  surfaces  at  each  time  step  and  are  useful  for  assessing 
iterative  convergence  of  integral  variables.  Forees  and  moments  are  broken  down  into  9 
eomponents  with  eontributions  due  to  skin  frietion,  piezometrie  pressure,  and  hydrostatie 
pressure  in  eaeh  of  the  3  eoordinate  direetions.  FNAMEO.  con v  and  FNAMEO .  f  scon v 
eontain  iterative  eonvergenee  history  of  the  point  variables  {U,  P,  Ur)  and  free-surfaee  wave 
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elevation  Q  respectively.  Writing  frequency  is  specified  by  the  namelist  variable 
it_save_conv,  which  is  otherwise  set  to  a  default  value  of  500  time  steps  (or  global 
iterations). 

The  remaining  output  file  is  the  global  solution  file,  which  contains  all  independent  and 
dependent  flow  variables.  By  default,  the  solution  is  written  only  on  the  last  time  step  for  steady 
flows  or  every  500  time  steps  for  unsteady  flows.  If  solutions  are  required  more  frequently,  for 
example  to  construct  animations,  the  namelist  variable  it_save_tec  can  be  used  to  specify 
desired  time-step  frequency.  This  file  is  formatted  as  an  ASCII  Tecplot  file  and  therefore  is 
directly  readable  by  the  commercial  visualization  software  Tecplot. 


Table  11.  Output  files 


Filename 

Unit  Number 

Description 

FNAMEO.rsto 

16 

Restart  file 

ENAMEO.tec 

26 

Global  solution  file  in  TECPLOT  format 

ENAMEO.res 

36 

Solution  residuals 

ENAMEO.conv 

46 

Iterative  history  of  pressure,  friction  velocity,  and 
velocity  along  no-slip  surfaces 

ENAMEO. forces 

56 

Iterative  history  of  forces 

ENAMEO  .moments 

57 

Iterative  history  of  moments 

ENAMEO. fsconv 

66 

Iterative  history  of  free  surface 

ENAMEO  body  forces,  tec 

86 

Grid,  blanking  variable  indicating  whether  point  is 
in/out  of  propeller  disk,  and  body-force 
components. 

7.  RECOMMENDED  VERIFICATION  AND  VALIDATION  PROCEDURES 

CFD  is  fast  becoming  an  integral  tool  in  the  engineering  design  process  as  it  is  applied  to 
increasing  complex  geometry  and  physics.  As  with  use  of  experimental  fluid  dynamics  (EFD)  in 
making  design  decisions,  assessment  of  quality  of  results  is  imperative,  which  has  accelerated 
progress  on  development  of  verification  and  validation  (V&V)  methodology  and  procedures  for 
estimating  numerical  and  modeling  errors  and  uncertainties  in  CFD  simulations.  However,  in 
spite  of  the  progress,  the  various  viewpoints  have  not  yet  fully  converged  and  current 
methodology  and  procedures  are  not  yet  standardized. 

Here,  however,  the  recommended  V&V  procedures  are  those  provided  by  Stern  et  al. 
(2001)  and  for  which  Wilson  et  al.  (2001)  presented  a  detailed  case  study  for  a  RANS  simulation 
of  an  established  benchmark  for  ship  hydrodynamics.  The  methodology  and  procedures 
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presented  therein  provides  a  pragmatic  approach  for  estimating  simulation  errors  and 
uncertainties.  The  philosophy  is  strongly  influenced  by  EFD  uncertainty  analysis.  The  approach 
allows  for  treatment  of  simulation  errors  as  either  stochastic  or  deterministic  and  properly  takes 
into  account  uncertainties  in  both  the  simulation  and  the  data  in  assessing  the  level  of  validation. 
A  brief  summary  of  the  methodology  and  procedures  is  provided  in  the  following. 


7.1  Methodology 

The  simulation  error  5^  is  defined  as  the  difference  between  a  simulation  result  S  and  the 
truth  T  and  is  composed  of  modeling  and  numerical  errors  {5g  =S-T  =  dgj^  +<^5Ar) 

with  corresponding  simulation  uncertainty  given  by  =  •  For  certain  conditions, 

both  the  sign  and  magnitude  of  the  numerical  error  can  be  estimated  as  (5’^^  =  ^ where 

dlf^  is  an  estimate  of  the  sign  and  magnitude  of  (5’^^  and  Ssn  is  the  error  in  that  estimate.  The 
simulation  value  is  corrected  to  provide  a  numerical  benchmark  Sc,  which  is  defined  by 

5, =5-^;^  (89) 

with  error  equation  -T  -  3^,^  +  and  corresponding  uncertainty  equation 

^lc~  ^SM  +  ^ScN  where  is  the  uncertainty  in  the  corrected  simulation  and  is  the 

uncertainty  estimate  for  ssn- 

Verification  is  defined  as  a  process  for  assessing  simulation  numerical  uncertainty  and, 
when  conditions  permit,  estimating  the  sign  and  magnitude  of  the  simulation  numerical 
error  itself  and  the  uncertainty  in  that  error  estimate  Ug^g; .  Numerical  error  is  decomposed  into 
contributions  from  iteration  number  Sj ,  grid  size  ,  time  step  dj. ,  and  other  parameters  dp , 
which  gives  the  following  expression  for  the  simulation  numerical  uncertainty 

Ugg,  =  uf  +Ul+Up+Ul  (90) 

For  situations  when  the  solution  is  corrected  to  produce  a  numerical  benchmark  the 
estimated  simulation  numerical  error  (5’^^  and  corrected  uncertainty  17 5^^  are  given  by 


dgjg  —  dj  +  Sq  +  dj  +  dp 


=u^  +u^  +u^  +u^ 

^  ScN  Gc  ^  Tc  ^  Pc 
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(91) 

(92) 


Validation  is  defined  as  a  proeess  for  assessing  simulation  modeling  uneertainty  by 
using  benehmark  experimental  data  and,  when  eonditions  permit,  estimating  the  sign  and 
magnitude  of  the  modeling  error  5^,^^  itself.  The  eomparison  error  E  is  given  by  the  differenee  in 
the  data  D  and  simulation  S  values 

E  =  D-S  =  S,-{S,^,+S,,^+S,,)  (93 ) 

where  has  been  deeomposed  into  the  sum  of  Sspd,  error  from  the  use  of  previous  data  such 

as  fluid  properties,  and  Ssma,  error  from  modeling  assumptions.  To  determine  if  validation  has 
been  achieved,  E  is  compared  to  the  validation  uncertainty  Uv  given  by 

Ul-Ul+Ul,.UU  (94) 

If  \E\  <  Uy ,  the  combination  of  all  the  errors  in  D  and  S  is  smaller  than  Uv  and  validation  is 
achieved  at  the  Uy  level.  If  Uv«\E\,  the  sign  and  magnitude  of  E=Ssma  can  be  used  to  make 

modeling  improvements.  For  the  corrected  approach,  the  equations  equivalent  to  equations  (93) 
and  (94)  are 

Ec  =  D  —  Sc  =  Sp—  +  dgpc  +  )  (95) 

K  =  K  -  ^Ima =ul+ f/L + Kn  (96) 


7.2  Procedures 

The  overall  CFD  V&V  procedures  can  be  conveniently  grouped  into  four  consecutive 
steps:  preparation,  verification,  validation,  and  documentation. 

Verification  is  accomplished  through  parameter  convergence  studies  using  multiple  solutions 
(at  least  3)  with  systematic  parameter  refinement  by  varying  the  input  parameter  while 
holding  all  other  parameters  constant.  Iterative  errors  must  be  accurately  estimated  or  negligible 
in  comparison  to  errors  due  to  input  parameters  before  accurate  convergence  studies  can  be 
conducted.  Changes  between  medium-fine  and  coarse-medium  ~  ~  ^k2 

solutions  are  used  to  define  the  convergence  ratio 

Rk  =  ^21  (97) 
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and  to  determine  convergenee  condition  where  ,  5*^^ ,  5*^^  correspond  to  solutions  with  fine, 

medium,  and  coarse  input  parameter,  respectively,  corrected  for  iterative  errors.  Three 
convergence  conditions  are  possible: 


(/)  Monotonic  convergence:  0  <  <1 

(n)  Oscillatory  convergence:  <  0  (98) 

(in)  Divergence: 

For  condition  (/),  generalized  RE  is  used  to  estimate  17^  or  (7*  and  17 .  For  condition  (n), 
uncertainties  are  estimated  simply  by  attempting  to  bound  the  error  based  on  oscillation 


maximums  Su  and  minimums  Sl,  i.e.. 


For  condition  (iii),  errors  and 


uncertainties  cannot  be  estimated. 

For  convergence  condition  (/),  generalized  RE  is  used  to  estimate  the  error  due  to 
selection  of  the  kth  input  parameter  and  order-of-accuracy  pt 


^RE. 


'2h 


Pk 


r/‘-l 


(99) 


(100) 


In(r^) 

Correction  of  equation  (99)  through  a  multiplication  factor  Ck  accounts  for  effects  of  higher- 
order  terms  and  provides  a  quantitative  metric  to  determine  proximity  of  the  solutions  to  the 
asymptotic  range 


o  *  o  * 


'2h 


-1 


(101) 


where  the  correction  factor  is  given  by 

-1 

C,=— -  (102) 

-1 

and  Pi^  is  an  estimate  for  the  limiting  order  of  accuracy  as  spacing  size  goes  to  zero  and  the 
asymptotic  range  is  reached  so  that  ^  1 .  When  solutions  are  far  from  the  asymptotic  range. 
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Q  is  sufficiently  less  than  or  greater  than  1 ,  only  the  magnitude  of  the  error  is  estimated  through 


the  uncertainty  U 

u,=\c,d''^]+\{\-c,)d;,]  (103) 

When  solutions  are  close  to  the  asymptotic  range,  is  close  to  1  so  that  51  is  estimated  using 
equation  (101)  and  Uf.  is  estimated  by 


= 


(104) 


Alternatively,  a  factor  of  safety  approach  proposed  in  Roach  (1998)  can  be  used  to  define 


and  U.  . 

Validation  is  accomplished  through  comparisons  with  benchmark  EFD  data,  including 
experimental  uncertainty  estimates  Ud-  If  the  three  variables  Uv,  \E\,  and  Ureqd  (programmatic 
validation  requirement)  are  considered,  there  are  six  combinations.  For  three  cases,  \E\<Uv  and 
validation  is  achieved  at  the  Uv  level,  but  for  only  one  of  these  Uv  <Ureqd  so  that  validation  is 
also  achieved  at  Ureqd-  In  these  cases,  attempting  to  estimate  modeling  errors  5sma  is  not  feasible 
from  an  uncertainty  standpoint.  For  the  three  other  cases,  Uv  <|F'|  and  using  the  sign  and 
magnitude  of  E  to  estimate  5sma  is  feasible  from  an  uncertainty  standpoint.  In  one  of  these  cases, 
Uv<\E\<Ureqd  SO  that  Validation  is  successful  at  the  \E\  level  from  a  programmatic  standpoint. 
Similar  conclusions  can  be  reached  using  the  corrected  comparison  error  and  corrected  validation 
uncertainty. 


8.  EXAMPLE  SIMULATION:  OPEN-WATER  PROPELLER  P5168 

In  this  section,  an  example  simulation  for  DTMB  open- water  propeller  P5168  is 
presented  and  discussed.  The  intention  is  to  demonstrate  some  of  the  capabilities  of  the  code 
including  overset  gridding,  non-inertial  relative  frames,  and  detailed  high-fidelity  resolution  of 
both  geometry  and  physics.  Discussion  follows  the  CFD  process  defined  in  Figure  1.  Input  files 
for  both  the  example  and  other  training  problems,  including  some  with  free-surface  effects,  may 
be  downloaded  from  the  CFDSHIP-IOWA  website,  http://www.iihr.uiowa.edu/~cfdship. 


62 


8.1  Geometry,  Benchmark  Data,  and  Conditions 

DTMB  propeller  model  P5168,  shown  in  Figure  22,  is  a  five-bladed,  controllable  pitch 
propeller  with  a  design  advance  ratio  of  J=\.21.  Chesnakas  and  Jessup  (2000)  presented  detailed 
velocity  field  LDV  measurements,  which  were  made  in  the  NSWC-CD  36-mch  water  tunnel. 
Their  study  was  undertaken  as  part  of  a  joint  project  with  the  Royal  Netherlands  Navy  and 
Marine  Research  Institute  to  develop  propeller  blade  tip  geometries  that  display  improved  tip- 
vortex  cavitation  characteristics.  Due  to  the  well-documented  experiment,  which  includes  data 
uncertainties,  P5168  has  become  an  international  benchmark  that  has  been  extensively  used  for 
CFD  validation.  Relevant  example  is  the  work  of  Chen  (2000)  and  Hsiao  and  Pauley  (1999), 
both  of  whom  used  P5168  to  undertake  detailed  study  of  propeller  tip  vortices  and  impact  of  grid 
resolution  and  non-linear  turbulence  closures.  In  this  report,  P5I68  is  used  to  demonstrate 
capability  of  CFDSHIP-IOWA  v3.03  in  simulating  propulsor  hydrodynamics  including  relative- 
frame  formulation,  near-wall  turbulence  models,  and  chimera-overset  gridding. 


(a)  View  looking  upstream  (b)  Side  view 

Figure  22.  P5168  Geometry. 

P5168  geometry  is  shown  in  Figure  22.  It  was  obtained  from  NSWC-CD  in  the  form  of 
an  IGES  file.  The  simulated  geometry  uses  an  infinite  shaft  of  constant  radius  whereas  the 
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experimental  geometry  ineluded  a  eylindrieal  fairwater  that  extended  96.8mm  downstream  of  the 
hub.  In  addition,  the  tested  geometry,  as  shown  in  Figure  1  of  Chesnakas  and  Jessup  (2000),  had 
a  small  inerease  in  hub  radius  near  the  propeller,  a  feature  whieh  is  not  ineluded  in  the  CFD 
model.  Sinee  the  foeus  here  is  simulation  at  near-design  operating  eonditions,  it  is  assumed  that 
the  inflow  is  uniform  in  the  eireumferential  direetion  and  that  the  flow  is  steady.  This  allows  the 
eomputational  domain  to  be  redueed  to  a  single  blade  passage  through  the  use  of  rotationally 
periodie  boundary  eonditions.  The  overall  domain,  whieh  was  shown  in  Figure  2,  extends  0.65D 
upstream,  1 .5D  downstream,  and  1 .5D  outward  in  the  radial  direetion. 

The  flow  eonditions  and  fluid  properties  are  based  upon  the  experimental  eonditions  at 
J=l.l,  whieh  is  slightly  off  the  design  point  and  is  used  sinee  it  results  in  a  stronger  tip  vortex 
than  J=1.27.  This  eondition  also  eorresponds  to  n=1450  rpm,  f/=10.70m/s,  7'=4715N,  g=481N- 
m.  Non-dimensional  parameters  used  in  CFDSHIP-IOWA  are  7?e=3.0xl0^  and  Ox=-27r/J=- 
5.712.  Flow  eonditions  are  speeified  in  the  namelist  input  file,  cfd_ship .  nml.  The  file 
used  for  P5168  is  ineluded  in  the  Appendix  B.l  as  an  example. 

For  detailed  propulsor  simulations  (as  opposed  to  using  a  propeller  body  foree),  the  only 
model  ehoiee  is  that  of  turbulenee  model.  For  the  simulations  here,  the  standard  blended  k-co/k-s 
turbulenee  model  was  used,  i.e.,  itm=2,  itm_switch=0. 

Steady  flow  simulations  use  free-stream  initial  eonditions  where  all  flow  variables  are  set 
to  UINF=1  and  VINF=WINF=0.  This  is  aehieved  by  setting  mode=0  in  the  cfd_ship .  nml 
fide.  As  shown  in  Figure  2,  the  following  boundary  eonditions  were  used  for  a  single-blade 
propeller  simulation:  relative-frame  no-slip  (ibtyp=22)  on  the  blade  and  hub  surfaees,  preseribed 
veloeity  (ibtyp=14)  on  the  upstream  inlet  plane,  exit  (ibtyp=12)  on  the  downstream  exit  plane, 
far-field  (ibtyp=13)  on  the  outer  boundaries,  and  rotationally-periodie  (ibtyp=52)  on  the  periodie 
faees.  The  preseribed  veloeity  profile  is 


U{r) 


(105) 


where  r/,  =  0.122199  is  the  hub  radius  and  r<y=  0.185  is  the  boundary  layer  thiekness  at  the  inlet 
plane.  Inlet  profile  in  equation  (105)  was  determined  by  Chen  (2000)  to  give  proper  loading  at 
the  blade  root.  Remaining  bloek  boundaries  are  either  patehed  multi-bloek  (ibtyp=92)  or 
Chimera  outer  boundaries,  the  latter  of  whieh  requires  use  of  PEGASUS  v5.1  from  NASA  Ames 
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and  which  will  be  further  diseussed  in  the  next  paragraph.  An  example  boundary  eondition  file 
is  ineluded  in  the  Appendix  B.3. 

8.2  Computational  Grids  and  Input  Parameters 

Generating  high-quality  patched-multi-block  grids  for  open-water  propellers  is 
ehallenging  due  radial  piteh  distribution,  espeeially  beyond  the  tip,  requirements  for  rotational 
periodieity,  and  need  for  boundary-layer  resolution  near  blade  and  hub  surfaees.  Together,  these 
issues  make  it  diffieult  to  eontrol  both  skewness  and  expansion  ratios.  CFDSHIP-IOWA  v3.03 
is  fairly  sensitive  to  grid  quality  sueh  that  grids  used  in  earlier  versions  (espeeially  those  based 
upon  the  Finite  Analytic  method)  often  result  in  unstable  simulations. 


(a)  Passage  bloeks  (b)  Blade  bloeks 

Figure  23.  P5168  Overset  Grid  System. 


Although  a  best  praetiees  document  for  CFDSHIP-IOWA  does  not  yet  exist  (e.g.,  a 
doeument  similar  to  Chan  et  ah,  2002),  the  following  guidelines  can  be  provided.  As  a  rule  of 
thumb,  CFDSHIP-IOWA  gives  most  accurate  results  using  expansion  ratios  less  than  1.5  and,  in 
general,  values  greater  than  3  should  be  avoided  sinee  simulation  stability  is  sensitive  to  this 
parameter.  At  multi-bloek  interfaees,  spaeing  on  eaeh  side  of  the  interfaee  should  be  similar. 
For  boundary-layer  resolution,  near-wall  spaeing  should  be  set  sueh  that  the  non-dimensional 
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wall  coordinates^  =  u^yjv  is  less  than  2.5  and  ideally  near  1.0.  Overset  gridding  provides  an 
approach  to  achieve  these  metrics  and  obtain  high-quality  grids,  especially  for  complex 
geometries. 

The  overset  grid  system  used  for  P5168  is  shown  in  Figure  23.  This  system  is  based 
upon  a  series  of  H-type  “passage”  blocks,  which  are  generated  without  regard  for  the  blade 
geometry,  and  0-type  “blade”  blocks,  which  are  generated  without  regard  for  the  periodic 
surfaces. 


Table  11.  P5 168  Grid  Parameters 


fine 

medium 

coarse 

Block  #  Name 

imax  jmax  kmax 

total 

imax  jmax  kmax 

total 

imax  jmax  kmax 

total 

1 

passage 1 

61 

81 

81 

400,221 

43 

58 

58 

144,652 

31 

41 

41 

52,111 

2 

passage2 

61 

81 

81 

400,221 

43 

58 

58 

144,652 

31 

41 

41 

52,111 

3 

passage3 

61 

81 

81 

400,221 

43 

58 

58 

144,652 

31 

41 

41 

52,111 

4 

passage4 

61 

81 

81 

400,221 

43 

58 

58 

144,652 

31 

41 

41 

52,111 

5 

inner_p 

41 

61 

31 

77,531 

29 

43 

22 

27,434 

21 

31 

16 

10,416 

6 

inner  s 

41 

61 

31 

77,531 

29 

43 

22 

27,434 

21 

31 

16 

10,416 

7 

tip_pl 

45 

61 

31 

85,095 

32 

43 

22 

30,272 

23 

31 

16 

11,408 

8 

tips  1 

45 

61 

31 

85,095 

32 

43 

22 

30,272 

23 

31 

16 

11,408 

9 

tip_p2 

45 

61 

31 

85,095 

32 

43 

22 

30,272 

23 

31 

16 

11,408 

10 

tip_s2 

45 

61 

31 

85,095 

32 

43 

22 

30,272 

23 

31 

16 

11,408 

11 

tip_p3 

41 

61 

45 

112,545 

29 

43 

32 

39,904 

21 

31 

23 

14,973 

12 

tip_s3 

41 

61 

45 

112,545 

29 

43 

32 

39,904 

21 

31 

23 

14,973 

13 

root_pl 

45 

61 

65 

178,425 

32 

43 

46 

63,296 

23 

31 

33 

23,529 

14 

root_p2 

41 

61 

65 

162,565 

29 

43 

46 

57,362 

21 

31 

33 

21,483 

15 

root_p3 

45 

61 

65 

178,425 

32 

43 

46 

63,296 

23 

31 

33 

23,529 

16 

root  si 

45 

61 

65 

178,425 

32 

43 

46 

63,296 

23 

31 

33 

23,529 

17 

root  s2 

41 

61 

65 

162,565 

29 

43 

46 

57,362 

21 

31 

33 

21,483 

18 

root_s3 

45 

61 

65 

178,425 

32 

43 

46 

63,296 

23 

31 

33 

23,529 

19 

phantom 

121 

51 

11 

67,881 

121 

51 

11 

67,881 

121 

51 

11 

67,881 

Total 

3,360,246 

Total 

1,202,280 

Total 

441,936 
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(a)  coarse 

Figure  24. 


(b)  medium  (c)  fine 

Systematic  grid  refinement,  surface  meshes. 


(a)  coarse  (b)  medium  (c)  fine 

Figure  25.  Systematic  grid  refinement,  cut  at  x/D=0.07. 


(a)  coarse 


(b)  medium 


(c)  fine 


Figure  26.  Systematic  grid  refinement,  overset-grid  system  at  trailing-edge  root. 
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This  approach  improves  grid  quality,  as  measured  by  orthogonality  and  expansion  ratios, 
aeeelerates  grid-generation  proeess,  and  permits  the  use  of  overset  refinement  meshes  (e.g.,  Kim 
et  ah,  2003).  The  grid  system  is  generated  using  GRIDGEN  from  Pointwise,  Ine.  and  was 
designed  for  a  3-grid  verifieation  study  as  shown  in  Table  1 1  and  with  near-wall  spaeing  less 
than  5x10'^  on  all  grids.  The  fine  grid  was  generated  using  GRIDGEN.  The  medium  and  eoarse 
grids  were  obtained  from  the  fine  grid  using  a  V2  refinement  ratio  in  eaeh  eoordinate  direetion. 
A  grid-sequeneing  tool  based  upon  linear  interpolation  in  eomputational  spaee  is  used.  This 
results  in  a  eoarse  grid,  whieh  is  systematieally  similar  to  the  fine  grid  due  to  the  fact  it  is  created 
through  grid  halving.  In  eontrast,  the  medium  grid  is  not  exaetly  systematieally  similar  to  the 
fine  grid  due  to  handwork  in  GRIDGEN  required  to  eorreet  geometry  errors  introdueed  in  using 
linear  interpolation.  Eigure  24  shows  the  surfaee  grids  for  each  grid  system.  Total  number  of 
grid  points  range  from  3.36  to  0.44  million  points  for  the  fine  and  eoarse  grids,  respectively. 

Before  running  CEDSHIP-IOWA,  overset-grid  interpolation  eoefficients  must  be 
eomputed  using  overset-grid  eommunieation  software.  Execution  of  PEGASUS  version  5.1  is 
aeeomplished  following  the  instruetions  in  Seetion  4.7  of  this  report  and  using  the  sample  input 
file  found  in  Appendix  B.2.  Therein,  it  ean  be  seen  that  a  novel  approaeh  has  been  taken  where 
two  separate  hole  eutters  ($HCUT  NAME)  have  been  defined.  The  first  one  is  a  traditional 
external-type  HCUT  named  “blade_cutter.”  This  HCUT  ereates  a  hole  in  the  passage  bloeks  due 
to  the  blade  surfaees  and  uses  a  phantom  mesh  at  the  root  of  the  blade  so  as  to  ereate  the  required 
“leak-free”  surfaee.  The  seeond  one  is  an  internal-type  hole  cutter  named  “passage  cutter” 
whose  purpose  is  to  trim  away  points  in  the  blade  bloeks  that  extend  past  the  periodie  surfaees. 
Eigures  25  and  26  show  the  eomposite  grid  system  for  coarse,  medium,  and  fine  grids.  Note  that 
the  fringe  boundaries  are  a  funetion  of  grid  resolution.  This  is  due  to  Eevel-2  interpolation. 
Impaet  on  verifieation  grid  studies  is  unknown. 

Numerieal  parameters  are  speeified  in  the  cfd_ship .  nml  file  of  whieh  an  example  for 
P5168  ean  be  found  in  the  Appendix  B.l.  Therein,  it  ean  be  seen  that  the  following  numerieal 
parameters  were  speeified:  2“‘*-order  upwind  spatial  aeeuraey  (ispat_order=3);  U*-order 
baekward  temporal  aeeuraey,  or  steady  flow  (itemp_order=l);  Cartesian  relative-frame 
coordinate  system  (icoord=2);  time  step  of  delt=0.01,  free  stream  veloeities  eomponents  of 
uinf=l  and  vinf=winf=0;  starting  and  ending  iterations  of  its=l  and  itend=10000; 
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iterative  frequency  of  writing  convergence  history  (it_save_conv),  restart  file 
(it_save_rst),  and  tecplot  file  (it_save_tec)  of  50,  500,  and  5000,  respectively;  number 
of  momentum  (ituvw),  velocity-pressure  coupling  (itvpc),  pressure  (itpr),  and  turbulence 
(itturb)  sub-iterations  of  5,  3,  5,  and  5,  respectively;  relaxation  factors  of  rfv=0.2, 
rfvb=0.2,  rfp=0.1,  and  rfpb=0.05;  fully  half-cell  formulation  of  the  pressure  equation 
(gama_pr=1.0);  line-solvers  are  turned  on  for  each  coordinate  direction  iswpl=l,  iswp2=l, 
iswp3=l;  pressure  reference  coordinate  location  is  specified  to  be  mref=l,  iref=l, 
j ref =41,  kref=21;  and  convective  discretization  of  the  k  and  co  equations  is  set  to  fi*-order 
upwind  (ispat_order_tm=l). 

8.3  Computing  Platforms 

Simulations  were  undertaken  using  two  computer  systems,  the  512-processor  SGI  Origin 
3800  at  the  Army  Research  Laboratory  Major-Shared  Resource  Center  (ARL-MSRC, 
http://www.arl.hpc.mil/)  and  the  72-processor  Linux  Beowulf  cluster  at  the  Penn  State  Applied 
Research  Laboratory.  The  former  is  a  distributed-shared  parallel  computer,  which  is  capable  of 
using  mixed-mode  parallelism  (i.e.,  code  is  compiled  as  make  02K_MP  I_OMP)  and  the  latter  is 
a  distributed  parallel  computer,  which  is  capable  of  MPI  parallelism  only  (i.e.,  code  is  compiled 
as  make  CLUSTER_MPI).  As  shown  in  Table  11,  block-size  distribution  is  not  uniform,  i.e., 
passage  blocks  are  a  factor  of  5  larger  than  the  smallest  block.  Therefore,  on  the  SGI,  some 
degree  of  load  balancing  was  achieved  through  the  use  of  OMP  threads,  which  is  an  automatic 
process  if  total_num_proc  is  specified  in  the  cfd_ship .  nml  file.  Setting 
total_num_proc  =  38,  4  OpenMP  threads  were  used  on  the  passage  blocks,  1  thread  on 
blocks  5-10,  and  2  threads  on  the  remaining  blocks.  The  computational  rate  was  5.5x10'^ 
processor-secs/grid-point/iteration.  For  comparison,  on  the  Linux  cluster,  where  OMP  threads 
cannot  be  used  for  load  balancing,  computational  rate  was  6.9x10“'^  processor-secs/grid- 
point/iteration. 

8.4  Verification  and  Validation  Results 

Typically,  iterative  convergence  is  assessed  using  both  force  and  moment  histories  and 
residuals  based  upon  the  change  in  variable  between  iterations.  However,  in  simulations  with 
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overlapping  surface  grids,  iterative  history  of  forces  and  moments  are  not  currently  available  due 
to  the  lack  of  a  run-time  interface  with  FOMOCO.  As  discussed  in  Section  4.8,  FOMOCO, 
which  is  a  part  of  the  suite  of  Chimera  Grid  Tools,  computes  forces  and  moments  on  overlapping 
grids  only  as  a  post-processing  step.  Therefore,  iterative  convergence  is  demonstrated  through 
the  use  of  residuals.  Coarse  grid  simulation  was  run  for  15,000  iterations  and  shows  4  orders 
magnitude  drop  for  all  variables.  A  grid  sequencing  approach  was  used  wherein  coarse-  and 
medium-grid  solutions  were  used  as  medium-  and  fine-grid  simulation  initial  conditions, 
respectively.  Medium  simulation  was  run  5000  iterations  and  fine  simulation  was  run  1000 
iterations.  Based  upon  previous  experience,  it  is  assumed  that  Ui  is  approximately  zero  for  all 
grids. 

Table  12.  P5168  Thrust  and  Torque  Coefficients  and  Comparison  to  Data. 


EFD  Data 

Coarse 

Medium 

Fine 

Kj 

0.313 

0.310 

0.318 

0.322 

E%=(D-S)/DM00 

1.0% 

-1.6% 

-2.9% 

s 

0.008 

0.004 

IOKq 

0.783 

0.765 

0.793 

0.805 

E%=(D-S)/DM00 

2.3% 

-1.3% 

-2.8% 

s 

0.028 

0.012 

Comparing  thrust  and  torque  coefficients  to  EFD  data  and  calculating  changes  between 
grids  s,  as  done  in  Table  12,  it  is  shown  that  near-blade  integral  quantities  display  monotonic 
grid  convergence  as  indicated  by  the  convergence  ratios  of  Rg  =  0.5  and  Rg  =  0.42  for  thrust  and 
torque,  respectively.  Based  upon  the  rules  of  Equation  (99),  grid  uncertainty  can  be  estimated 
using  RE  and  is  shown  in  Table  13.  Since  Ud  is  not  available  for  Kj  and  Kq,  validation 
uncertainty  Uyi^s  not  calculated. 


Table  13.  Verification  of  P5168  Thrust  and  Torque  Coefficients. 


PG 

RE 

Cg 

Ugc 

Ui 

UsN 

Kj 

2.0 

0.004 

1.0 

0.0 

0.0 

0.004 

Kq 

2.4 

0.015 

1.29 

0.0043 

0.0 

0.015 
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(a)  coarse  (b)  medium 

Figure  27.  Surface  pressure  on  pressure  side  of  blade. 


(c)  fine 


(a)  coarse  (b)  medium  (c)  fine 

Figure  28.  Surface  pressure  on  suction  side  of  blade. 


Figure  29.  Comparison  of  surfaee  pressure  at  r/Rp=0.716. 
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Surface  pressure  contours  are  shown,  for  each  grid  system,  in  Figures  27  and  28.  In 
general,  typical  blade  pressure  distributions  are  shown  with  maximum  loading  at  approximately 
y4  span,  high  magnitude  on  the  pressure-side  leading-  and  trailing-edges,  and  low  magnitude  on 
the  suction-side  tip.  Contours  also  display  monotonic  grid  convergence  with  fine  grid  displaying 
highest  resolution  of  detail.  Figure  29  shows  a  comparison  of  surface  pressure  at  r/Rp=0.716. 
Pressure-side  of  blade  displays  monotonic  grid  convergence  whereas  suction  side  displays  more 
oscillatory  behavior  along  the  chord. 

Figure  30  highlights  the  grid  convergenee  of  the  tip  and  root  vortices  using  iso-surfaees 
of  intrinsic  swirl  parameter  (Berdahl  and  Thompson,  1993)  at  level  t  =  8.0  colored  by 
normalized  helicity  H  =  coxU ,  the  latter  of  which  indicates  direction  of  rotation.  Counter¬ 
rotating  root  vortices  display  inereasing  strength,  as  indieated  by  inereased  persistence  and 
organization  in  the  downstream  direction,  with  increasing  grid  resolution.  Similar  observation 
can  be  made  for  tip  vortex,  i.e.,  fine  grid  shows  tip  vortex  with  longest  persistenee,  however,  the 
very  short  persistence  of  the  coarse  grid  suggests  grid  in  tip-vortex  region  is  not  in  the 
asymptotic  range. 


(a)  coarse  (b)  medium  (c)  fine 


Figure  30.  Tip-  and  root-vortex  visualization  using  iso-surface  of  intrinsic 
swirl  parameter  (t  =  8.0)  colored  by  normalized  helicity. 

Figure  31  shows  a  comparison  of  simulation  to  data  for  axial-velocity  contours  at 
x/D=0.1193.  In  general,  all  three  solutions  show  resolution  of  the  global  trends,  i.e.,  thin  blade 
wakes  with  increasing  wake  thickness  near  hub,  maximum  velocity  on  the  suction  side,  and 
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contours  typical  of  a  tip  vortex  at  r/D=0.5.  Moreover,  resolution  of  blade  wakes  and  tip-vortex 
show  grid  convergence  in  that  a  general  trend  of  improved  resolution  with  grid  can  be  observed. 
However,  detailed  comparison  of  certain  contour  levels  between  the  3  solutions  indicates  some 
degree  of  oscillatory  grid  convergence.  This  is  more  clearly  shown  in  Figure  32  which  shows 
the  axial  velocity  as  a  function  of  circumferential  position  at  a  constant  radius  of  r/D  =  0.465. 
This  figure  shows  the  blade-to-blade  variability  of  the  data.  Overall,  the  3  solutions  display  a 
non-monotonic  divergent  condition,  which  may  be  due  to  both  a  coarse  grid  solution  outside  of 
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the  asymptotic  range  and  grids  which  are  not  exactly  systematically  similar;  observations 
supported  by  Figure  30  and  previously  discussed  in  section  8.2,  respectively.  However,  medium 
and  fine  grids  appear  to  be  converging  towards  the  data  which  suggests  that  a  finer  grid  4th 
solution  at  approximately  9.5  million  grid  points,  which  would  be  a  very  large  grid  for  a  single¬ 
blade  simulation,  would  provide  3  solutions  in  the  asymptotic  range.  This  represents  the 
principal  challenge  of  RE  based  error  estimation  and  points  to  the  need  for  continued  research 
and  development  in  automatic  generation  of  geometrically  similar  overset  grids  and  interpolation 
coefficients,  and  single-grid  error  estimation  techniques  (e.g.,  Celik  et  ah,  2003). 


Figure  32.  Comparison  of  circumferential  distribution  of  axial  velocity  at  r/D=0.465  and 

x/D=0.1193. 
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9.  CONCLUDING  REMARKS 

CFDSHIP-IOWA  is  a  general-purpose  unsteady  Reynolds-averaged  Navier-Stokes  CFD 
eode  that  has  been  developed  to  handle  a  broad  range  of  ship  hydrodynamies  problems.  Purpose 
of  this  report  was  to  provide;  detailed  doeumentation  of  the  modeling,  numerical  methods,  and 
code  development;  user  instructions  on  creating  input  files  and  post-processing;  recommended 
verification  and  validation  procedures;  and  an  example  simulation.  As  a  framework  for 
achieving  successful  simulations,  an  approach  based  upon  formulation  of  an  initial  boundary 
value  problem  and  execution  of  a  well-defined  CFD  process  was  developed  and  followed 
throughout  the  report.  An  example  simulation,  and  other  recent  applications,  demonstrates  the 
capability  of  CFDSHIP-IOWA  v3.03  to  simulate  practical  ship  hydrodynamics  problems. 
Successful  use  in  both  thesis  and  project  research  and  transition  to  other  organizations 
demonstrates  the  success  of  the  overall  design  objectives. 

Largely  due  to  successes  such  as  those  shown  and  referred  to  herein,  role  of  CFD  in 
analysis  and  design  of  future  marine  vehicles  continues  to  expand.  Beyond  higher-fidelity 
simulation  of  resistance  and  propulsion  problems,  albeit  including  new  applications  and  off- 
design  conditions,  it  is  anticipated  that  future  simulations  will  move  towards  including 
environmental  effects  (e.g.,  seaway,  stratification,  shallow  water),  simulation-based  design  and 
optimization,  resolving  complex  maneuvering  scenarios  of  multiple-body  configurations  (e.g., 
submarine  or  surface-ship  launch  of  adjunct  vehicles),  and  supporting  improved  modeling  of 
acoustic  and  non- acoustic  signatures.  However,  in  spite  of  continued  advancements  in  HPC 
hardware,  the  magnitude  of  these  problems  are  expected  to  be  on  the  order  of  20-50  million  grid 
points  and  therefore  require  continued  development  of  more  accurate  numerics  and  faster 
computing  algorithms.  Areas  of  future  development  include  improved  algebraic  solvers, 
adaptive  gridding,  ideally  using  a  single-grid  error-estimation  equation,  implementation  of 
geometry  manipulation  protocols  for  specification  of  vehicle  configurations  and 
maneuvering/seakeeping  scenarios,  multi-phase  level-set  methods  for  robust  free-surface 
simulations,  and  hybrid  RANS-LES  models  for  numerous  applications. 
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APPENDICES 


APPENDIX  A:  FILE  FORMATS 

Detailed  format  deseription  of  input  and  output  fdes  are  deseribed  in  this  appendix. 
A.l  Grid  File 

The  grid  data  is  read  by  subroutine  get_grid  whieh  is  loeated  in  the 
cf  dship_stio  .  F  file.  The  file  format  is  as  follows. 


read(15,*)  nmesh 
do  m=l, nmesh 

read (15,*)  imax (m) ,  jmax(m),  kmax (m) 

enddo 

do  m=l, nmesh 

read(15,*)  ( ( (x ( I , j , k) , 1=1 , imax (m) ) , j =1 , jmax (m) , k=l , kmax (m) ) ,  & 

( ( (y (I, j , k) , 1=1, imax (m) ) , j=l, jmax (m) , k=l, kmax (m) ) ,  & 

(((z(I,j,k),I  =  l, imax (m) ) , j  =1 , jmax (m) , k=l , kmax (m) ) 

enddo 


A.2  Namelist  Input  File 

There  are  9  NAMELISTS  in  the  eode  and  each  must  appear  in  the  input  file, 
cf  d_ship  .  nml.  The  file  is  opened  and  read  in  subroutine  input_runtime  and 
subroutine  input_grid_variables ,  both  of  which  are  in  the  file 
cf dship_mods . F90 


open (unit=8 ,file='cfd  ship. nml ' , status= ' old ' , action= ' read ' , iostat=ierror ) 

read ( 8 , nml=control ) 

read ( 8 , nml=f low  parameters) 

read ( 8 , nml=grid  parameters) 

read ( 8 , nml=i deration) 

read ( 8 , nml=solver ) 

read (8, nml=turbulence) 

read ( 8 , nml=f ree_surface) 

read ( 8 , nml=propeller ) 

read ( 8 , nml=f ilenames ) 

close  ( 8 ) 


File  format  is  shown  in  Appendix  B.l,  however,  each  NAMELIST,  including  function  and 
variable  name,  description  and  default  values,  is  described  in  the  following. 
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$CONTROL 

This  input  sets  global  values 


Variable 

Description 

Default 

MODE 

Flag  to  set  initial  conditions  (=0  for  start  from  firee- 
stream,  =1  for  start  from  restart  file) 

0 

TOTAL  NUM  PROGS 

Total  number  of  processors  used  (only  used  for 
mixed-mode  parallelism) 

1 

I SPAT  ORDER 

Sets  spatial  order-of-accuracy  of  convective  terms 
(see  Table  4  for  options) 

3 

ITEMP  ORDER 

Sets  temporal  order  of  accuracy  (see  Table  3  for 
options) 

1 

I  COORD 

Coordinate  system  flag 

•  Cartesian,  absolute  frame,  ICOORD=l 

•  Cartesian,  relative  frame,  ICOORD=2 

•  Cylindrical,  absolute  frame,  ICOORD=3 

•  Cylindrical,  relative  frame,  ICOORD=4 

1 

$FLOW_PARAMETERS 
This  input  sets  flow  parameters. 


Variable 

Description 

Default 

RE 

Reynolds  number 

none 

ENUM 

Froude  number 

0.0 

DELT 

Time  step 

none 

TIME  RAMP  END 

Cubic  polynomial  ramp-up  time  of  unsteady  flow 

2.0 

UINE 

Free-stream  velocity  component  in  x-direction 

1.0 

VINE 

Free-stream  velocity  component  my  (or  r)  direction 

0.0 

WINE 

Free-stream  velocity  component  in  z  (or  0)  direction 

0.0 

$GRID_PARAMETERS 

This  input  sets  grid  modification  parameters.  All  variables,  except  for  the  ones  that  set  the  point 
about  which  moments  are  calculated,  are  arrays  that  can  have  a  unique  value  for  each  block  in 
the  grid  system. 


Variable 

Description 

Default 

X  TRANSLATE 

Distance  x-coordinate  is  translated.  Value  subtracted 
from  initial  grid 

0.0 

Y  TRANSLATE 

Distance  y-coordinate  is  translated.  Value  subtracted 
from  initial  grid 

0.0 

Z  TRANSLATE 

Distance  z-coordinate  is  translated.  Value  subtracted 
from  initial  grid 

0.0 
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ALPHA 

GAMA 

BETA 

SCALE_MESH 

X_ROT_CENT 

Y_ROT_CENT 

Z_ROT_CENT 

X_MOM_CENT 

Y_MOM_CENT 

Z_MOM_CENT 

AGVX 

AGVY 

AGVZ 


Rotation  about  z-axis  0.0 

Rotation  about  y-axis  0.0 

Rotation  about  x-axis  0.0 

Multiplicative  factor  to  scale  grid  1 .0 

x-coordinate  of  rotation  center  0.0 

y-coordinate  of  rotation  center  0.0 

z-coordinate  of  rotation  center  0.0 

x-coordinate  of  point  about  which  moments  calculated  0.0 
y-coordinate  of  point  about  which  moments  calculated  0.0 
z-coordinate  of  point  about  which  moments  calculated  0.0 
Angular  velocity  about  x-axis  0.0 

Angular  velocity  about  y-axis  0.0 

Angular  velocity  about  z-axis _ ^ 


Note  that  the  order  of  grid  manipulation  is  1)  scale,  2)  translate,  and  3)  rotate. 


$ ITERATION 

This  input  controls  iterative  solvers,  starting  and  ending  time  step,  and  convergence  tolerances. 


Variable 

Description 

Default 

ITS 

Starting  time  step  (or  global  iteration) 

1 

ITEND 

Ending  time  step  (or  global  iteration) 

None 

IT  SAVE  TEC 

Frequency  for  writing  tecplot  file 

500 

IT  SAVE  CONV 

Frequency  for  saving  convergence  history 

500 

IT  SAVE  RST 

Frequency  for  writing  restart  file 

500 

ITUVW 

Number  of  sub-iterations  for  solution  of  momentum 
equation 

5 

ITVPC 

Number  of  velocity-pressure  coupling  loops  (i.e.,  PISO 
pressure  correction  steps) 

2 

ITPR 

Number  of  sub-iterations  for  solution  of  pressure  equation 

5 

ITTURB 

Number  of  sub-iterations  for  solution  of  2-equation 
turbulence  model  equations 

5 

TOL  UVW 

Convergence  tolerance  for  momentum  equations  (used 
only  for  unsteady  flow) 

l.Oe-04 

TOL  PR 

Convergence  tolerance  for  pressure  equation  (used  only  for 
unsteady  flow) 

l.Oe-04 
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$ SOLVER 

This  input  sets  the  relaxation  parameters,  solver  sweep  direetions,  and  the  pressure-referenee 
loeation  for  Fr=0  simulations. 


Variable 

Description 

Default 

REV 

Velocity  relaxation  factor,  solver  level 

0.2 

REVB 

Velocity  relaxation  factor,  global  level  (=1.0  for 
unsteady  flows) 

0.2 

REP 

Pressure  relaxation  factor,  solver  level 

0.1 

REPB 

Pressure  relaxation  factor,  global  level  (=1.0  for 
unsteady  flows) 

0.1 

ISWPl, ISWP2, ISWP3 

Flag  to  set  active  directions  (^,ri,C  coordinates) 
for  line  solver 

0,1,0 

MREE, IREE, JREE, KREE 

Block  #  and  (i,j,k)  index  for  setting  reference 
pressure  location.  For  Fr=0,  pressure  at  this 
point  is  0.0 

1,1,  1,1 

$TURBULENCE 

This  input  selects  turbulence  model  and  options. 


Variable 

Description 

Default 

ITM 

Sets  turbulence  model. 

•  Laminar  flow,  ITM=0 

•  Baldwin-Lomax  model,  ITM=1 

•  Blended  k-®  model,  ITM=2 

0 

ITM  SWITCH 

Flag  to  set  model  options  for  Blended  k-co 
model 

•  Standard  model,  ITM_SWITCH=0 

•  SST  model,  ITM_SWITCH=1 

•  Wilcox  \ow-Re  model,  ITM_SWITCH=2 

0 

ITM  SPAT  ORDER  TM 

Sets  order-of-accuracy  of  2-eqn  turbulence 
model 

I SPAT  ORDER 
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$FREE_SURFACE 

This  input  sets  parameters  for  free-surfaee  solver  and  dynamie  grid  eonforming  proeess. 


Variable 

Description 

Default 

ITFSMAX 

Number  of  free-surfaee  iterations 

5 

TOL  FS 

Convergence  tolerance  for  free-surfaee  solution  (only  used  for 
unsteady  simulations) 

l.Oe-04 

WAVE LANK 

Blanking  distance  for  near-wall  region  of  free-surfaee 

0.0 

IFS 

Frequency,  in  global  time  steps,  of  free-surfaee  solution 

5 

ICFM 

Flag  indicating  whether  grid  is  conformed  to  free  surface  or 
design  waterline  (z/Z=0.0).  ICFM=0,  no  conform;  ICFM  =1, 
conform  turned  on. 

0 

$PROPELLER 

This  input  eontrols  the  preseribed  body-foree  deseribed  in  Seetion  2.5. 


Variable 

Description 

Default 

I  PROP 

number  of  propellers 

0 

CT 

Thrust  coefficient 

0.0 

CKQ 

Torque  coefficient 

0.0 

ADVANCE 

:  COEF 

Propeller  advance  coefficient 

0.0 

DXPROP 

Propeller  disk  thickness 

0.0 

RP 

Propeller  radius 

0.0 

X  PROP 

CENTER 

x-coordinate  of  propeller  center 

0.0 

Y  PROP 

CENTER 

y-coordinate  of  propeller  center 

0.0 

Z  PROP 

CENTER 

z-coordinate  of  propeller  center 

0.0 

RH 

Propeller  hub  radius  (in  decimal  %  of  RP) 

0.0 

SHAFTALPHA 

Angle  of  propeller  shaft  with  x-coordinate 

0.0 

$ El RENAMES 

This  input  sets  the  filename  extensions. 


Variable  Description _ Default 

EGRI D  Variable  which  sets  filename  for  grid  file  None 

FNAMEI  Variable  which  sets  “previous_simulation”  filename  prefix  None 

FNAMEO  Variable  which  sets  “current  simulation”  filename  prefix  None 
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A. 3  Boundary  Condition  File 


The  boundary  condition  data  is  read  by  subroutine  input_bcs  which  is  located  in 
the  cf  dship_stio  .  F  file.  The  detailed  file  format  is  as  follows. 

do  i=l,nmesh 

read  (iunit,*)  nbc(i) 
do  n=l , nbc  ( i ) 

read  (iunit,*)  ibtyp(n,i) 
read  (iunit,*)  ibdir(n,i) 
read  (iunit,*)  ibcs (n, i ) , ibce  (n, i ) 
read  (iunit,*)  jbcs (n, i) , jbce (n, i) 
read  (iunit,*)  kbcs (n, i ) , kbce (n, i ) 
read  (iunit,*)  ibcord(n,i) 

if (ibtyp (n, i) .eq. 30)  read  (iunit,*)  if sf ilter (n, i ) 
if ( ibtyp (n, i ) . eq . 91 . or . ibtyp (n, i ) .eq.92.or.  & 

ibtyp(n,i) . eq . 4 1 . or . ibtyp (n, i ) .eq.42.or.  & 

ibtyp(n,i) . eq . 51 . or . ibtyp (n, I ) .eq.52)  then 


read 

(iunit, *) 

ndmesh (n, i ) 

read 

(iunit, *) 

idbdir (n, i ) 

read 

(iunit, *) 

ides  (n,  i)  ,  idee  (n,  i 

read 

(iunit, *) 

jdes (n, i) , jdee (n,  i 

read 

(iunit, *) 

kdes  (n,  i)  ,  kdee  (n,  i 

endif 

enddo 

enddo 


A.4  Overset  Interpolation  Coefficient  File 

The  boundary  condition  data  is  read  by  subroutine  get_chimera  which  is  located 
in  the  cf  dship_chimera  .  F90  file.  The  detailed  file  format  is  as  follows. 


do  in=l,nmesh 

read (199)  ibpnts (m) , iipnts (m) , iieptr (m) , iisptr (m) ,  & 
imax  peg (m) , jmax  peg (m) , kmax  peg (m) 
read(199)  ( ii ( i ), i=iisptr (m) , iieptr (m) ) ,  & 

( j i  ( i ), i=iisptr (m) , iieptr (m) ) ,  & 

(ki (i)  ,  i=iisptr (m) , iieptr (m) ) ,  & 

(dxint_peg (i) , i=iisptr (m) , iieptr (m) ) ,  & 

(dyint_peg (i) , i=iisptr (m) , iieptr (m) )  ,  & 

(dzint_peg (i) , i=iisptr (m) , iieptr (m) ) 
read(199)  ( ib ( i ) , i=ibsptr (m) , ibeptr (m) ) ,  & 

( jb (i) , i=ibsptr (m) , ibeptr (m) ) ,  & 

(kb (i)  ,  i=ibsptr (m) , ibeptr (m)  ) ,  & 

(ibc (i) , i=ibsptr (m) , ibeptr (m) ) 
read (199)  (iblank_peg (i) , i=first (m) , last (m) ) 
enddo 
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A.J  Global  Solution  File 


The  global  solution  file  is  written  by  subroutine  save_tec  whieh  is  loeated  in  the 
cf  dship_stio  .  F  file.  The  detailed  file  format  is  as  follows. 


write (26, 10) 
write (26,20) 
do  m=l,nmesh 
write (26, 30) 
write (26, 40) 
write (26, 40) 
write (26, 40) 
write (26, 40) 
write (26, 40) 
write  (26,  40) 
write (26, 40) 
write (26, 40) 
write (26, 40) 
write (26, 40) 
write (26,41) 
write (26, 40) 
write (26, 40) 
enddo 


imax  (m)  ,  jmax  (m)  ,  kmax  (m) 

(  (  (x(i,j,k,m)  ,i=l,  imax  (m)  )  ,  j  =1 ,  jmax  (m)  )  ,  k=l ,  kmax  (m)  ) 

(  (  (y(i,j,k,m)  ,i=l,  imax  (m)  )  ,  j=l,  jmax  (m)  )  ,  k=l,  kmax  (m)  ) 
(((z(i,j,k,m),i=l,  imax  (m)  )  ,  j  =1 ,  jmax  (m)  )  ,  k=l ,  kmax  (m)  ) 

(  (  (u(i,j,k,m)  ,i=l,  imax  (m)  )  ,  j  =1 ,  jmax  (m)  )  ,  k=l ,  kmax  (m)  ) 

(  (  (v(i,j,k,m)  ,i=l,  imax  (m)  )  ,  j  =1 ,  jmax  (m)  )  ,  k=l ,  kmax  (m)  ) 

(  (  (w(i,j,k,m)  ,i=l,  imax  (m)  )  ,  j  =1 ,  jmax  (m)  )  ,  k=l ,  kmax  (m)  ) 

( ( (pr(i,j,k,m) ,i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l , kmax (m) ) 

( ( (zut(i,j,k,m) ,i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l , kmax (m) ) 

( ( (yplus (i, j , k,m) , i=l, imax (m) ) , j=l, jmax (m) ) , k=l, kmax (m) ) 

( ( (uplus(i,j,k,m) ,i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l , kmax (m) ) 

( ( (iblank(i,j,k,m) ,i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l , kmax (m) ) 
( ( (ak(i,j,k,m) ,i=l, imax (m) ) , j=l, jmax (m) ) , k=l, kmax (m) ) 

( ( (ao (i, j , k,m) , i=l, imax (m) ) , j=l, jmax (m) ) , k=l, kmax (m) ) 


10  format (' TITLE  =  "TITLE  STRING"') 

20  format ( 'VARIABLES  =  X,  Y,  Z,  U,  V,  W,  P,  ZUT,  YPLUS,  UPLUS,  IBLANK',', 
K,  OMEGA') 

30  format ('ZONE  I  =',I3,',  J=',I3,',  K=',I3,',  F=BLOCK') 

40  format ( 6el4 . 6) 

41  format (30i4) 


A.6  Restart  File 

The  restart  file  is  read  by  subroutine  get_restart  and  is  written  by 
subroutine  save_restart,  both  of  whieh  are  loeated  in  the  cf  dship_stio  .  F  file. 
The  detailed  file  format  is  as  follows. 


read (35)  ntot_orig,  nmesh 
do  m=l, nmesh 

read (35)  f ir st_orig (m) , 
read (35)  imax_orig (m) , 
enddo 

read(35)  (  (  (  (xO ( i, j , k, m) 
read (35) 
read  (35) 
read (35) 
read (35) 
read (35) 
read (35) 
read (35) 


( ( ( (yO (i, j , k,m) 

(  (  (  (zO  (i, j , k,m) , i 
(  (  (  (uO (i, j , k,m)  ,  i 
(  (  (  (vO (i, j , k,m) , i 
(  (  (  (wO (i, j , k,m) , i 
(  (  (  (prO (i, j , k,m)  , 
(  (  (  (zutO (i, j , k,m) 


last_orig (m) , 
jmax_orig (m) , 

=1, imax (m) ) , j= 
■1 ,  imax  (m)  )  ,  j  = 
■1 ,  imax  (m)  )  ,  j  = 
■1 ,  imax  (m)  )  ,  j  = 
■1 ,  imax  (m)  )  ,  j  = 
■1 ,  imax  (m)  )  ,  j  = 
i=l, imax (m) ) , j 
,  i=l ,  imax (m) ) , 


length_orig (m) 
kmax_orig (m) 

1, jmax (m) ) , k=l, 
1, jmax (m) ) , k=l, 
1,  jmax (m) ) , k=l, 
1, jmax (m) ) , k=l, 
1, jmax (m) ) , k=l, 
1,  jmax (m) ) , k=l, 
=1 , jmax (m) ) , k=l 
j=l, jmax (m) ) , k= 


kmax (m) ) , m=l , 
kmax (m) )  ,  m=l , 
kmax (m) ) , m=l , 
kmax (m) )  ,  m=l , 
kmax (m) ) , m=l , 
kmax (m) )  ,  m=l , 
, kmax (m) ) , m=l 
1 , kmax (m) ) , m= 


nmesh) 

nmesh) 

nmesh) 

nmesh) 

nmesh) 

nmesh) 

, nmesh) 

1 , nmesh) 
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read (35)  ( ( ( (akO (i, j , k,in) , i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((aoO(i,j,k,in),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

if(itemp  order. eq. 2)  then 

read (35)  ((((x00(i,j,k,m),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ( ( ( (yOO (i, j , k,m) , i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((z00(i,j,k,m),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((u00(i,j,k,m),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((v00(i,j,k,m),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((w00(i,j,k,m),i=l, imax (m) ) , j  =1 , jmax (m) ) , k=l, kmax (m) ) ,m=l , nmesh) 

read (35)  ((((akOO(i,j,k,m),i=l, imax (m) ) , j=l, jmax (m) ) , k=l, kmax (m) ) ,m=l, nmesh) 

read (35)  ((((ao00(i,j,k,m),i=l, imax (m) ) , j=l, jmax (m) ) , k=l, kmax (m) ) ,m=l, nmesh) 

endif 


APPENDIX  B:  SAMPLE  INPUT  FILES 

B.  1  CFDSHIP-IO  WA  Namelist  Input  File,  “cfd_ship.  nml  ” 

$control 
mode 

total_num_procs 
ispat_order 
itemp_order 
icoord 
$end 

$grid 
agvx 


$end 

$f low_parameters 
re 

delt 
uinf 
vinf 
winf 
$end 

$iteration 
its 
itend 

it_save_conv 
it_save_rst 
it_save_tec 
ituvw 
itvpc 
itpr 
itturb 
tol_uvw 
tol_pr 
$end 

$solver 
rfv  =  0.2 

rfvb  =  0.2 

rfp  =  0.10 

rfpb  =  0.05 

gama_pr  =  1.0 

iswpl  =  1 

iswp2  =  1 


=  00001 
=  10000 
050 
=  00500 
=  05000 
=  5 
=  3 
=  5 
=  5 

=  l.Oe-4 
=  l.Oe-4 


=  3.0e6 
=  0.01 
=  1.0 
=  0.0 
=  0.0 


parameters 

=  -5.712,  -5.712,  -5.712,  -5.712,  -5.712,  -5.712, 


-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5 . 712 

-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5.712 

-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5.712, 

-5.712 

=  0 
=  24 
=  3 
=  1 
=  2 
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iswp3  =  1 

mref  =  1 

iref  =  1 

jref  =  41 

kref  =  21 

$end 

$turbulence 
itm  =  2 

itm_switch  =  0 


ispat_order_tm  =  1 

$end 

$f ree_surf ace 
$end 

$propeller 

$end 

$f ilenames 

fgrid  =  "p5168_14b_cl . grd" 

fnamei  =  "p5168_14b_cl" 

fnameo  =  "p5168_14b_cl" 

$end 

B.  2  PEGASUS  v5. 1  Namelist  Input  File,  “peg.  in  ” 

$GLOBAL 

FRINGE  =  2, 

PROJECT  =.T., 

$END 


$MESH 

NAME 

= 

'BLK-1 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l,  $END 

$MESH 

NAME 

= 

'BLK-2 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-3 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-4 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-5 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-6 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-7 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-8 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-9 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-10 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-11 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-12 

r 

KINCLUDE= 

2, 

-1, 

OFFSET=l, $END 

$MESH 

NAME 

= 

'BLK-13 

r 

KINCLUDE= 

2, 

-1, 

LINCLUDE=  2, 

-1, 

$END 

$MESH 

NAME 

= 

'BLK-14 

r 

KINCLUDE= 

2, 

-1, 

LINCLUDE=  2, 

-1, 

$END 

$MESH 

NAME 

= 

'BLK-15 

r 

KINCLUDE= 

2, 

-1, 

LINCLUDE=  2, 

-1, 

$END 

$MESH 

NAME 

= 

'BLK-16 

r 

KINCLUDE= 

2, 

-1, 

LINCLUDE=  2, 

-1, 

$END 
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$MESH 

NAME  = 

'BLK-17  '  , 

KINCLUDE= 

2, 

-1, 

LINCLUDE= 

2,  -1, 

$END 

$MESH 

NAME  = 

'BLK-18  '  , 

KINCLUDE= 

2, 

-1, 

LINCLUDE= 

2,  -1, 

$END 

$MESH 

NAME  = 

'BLK-19 '  , 

KINCLUDE= 

2, 

-1, 

PHANTOM=.T, 

.  , $END 

$HCUT  NAME  =  ' blade_cutter ' , 

INTERNAL=.F.  , 

MEMBER  =  'BLK-19 ' , 

'BLK-5',  'BLK-6',  'BLK-7',  'BLK-8',  'BLK-9',  'BLK-10',  'BLK-11', 
'BLK-12',  'BLK-13',  'BLK-14',  'BLK-15',  'BLK-16',  'BLK-17',  'BLK-18', 
INCLUDE  =  ' BLK-1 '  ,  ' BLK-2 '  ,  ' BLK-3 '  ,  ' BLK-4 ' , 

CARTX=-81.6,  81.6, 

$END 

$HCUT  NAME  =  ' passage_cutter ' , 

INTERNAL=.T.  , 

MEMBER  =  ' BLK-1 '  ,  ' BLK-2 '  ,  ' BLK-3 '  ,  ' BLK-4 ' , 

INCLUDE  =  'BLK-5',  'BLK-6',  'BLK-7',  'BLK-8',  'BLK-9',  'BLK-10',  'BLK-11', 

'BLK-12',  'BLK-13',  'BLK-14',  'BLK-15',  'BLK-16',  'BLK-17',  'BLK-18', 

$END 

$LEVEL2  EXCLUDE= ' BLK-1 ' ,  ' BLK-2 '  ,  ' BLK-3 '  ,  ' BLK-4 '  , 

$END 


$BCINP  ISPARTOF  =  'BLK-1', 


IBTYP  = 

5, 

40, 

5, 

5, 

5, 

5 

IBDIR  = 

3, 

-1, 

-3, 

2, 

1, 

-2 

JBCS  = 

1, 

61, 

1, 

1, 

1, 

1 

JBCE  = 

t — 1 

61, 

t — 1 

t — 1 

1, 

61 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

81 

KBCE  = 

CO 

CO 

81, 

1, 

CO 

81 

LBCS  = 

1, 

1, 

81, 

1, 

1, 

1 

LBCE  = 
$END 

1, 

CO 

81, 

CO 

CO 

81 

$BCINP  ISPARTOF  = 

'BLK- 

2', 

IBTYP  = 

5, 

5, 

40, 

5, 

40, 

5 

IBDIR  = 

2, 

3, 

-1, 

-3, 

1, 

-2 

JBCS  = 

1, 

1, 

61, 

1, 

1, 

1 

JBCE  = 

t — 1 

61, 

61, 

t — 1 

1, 

61 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

81 

KBCE  = 

1, 

81, 

CO 

81, 

CO 

81 

LBCS  = 

1, 

1, 

1, 

81, 

1, 

1 

LBCE  = 
$END 

\ — 1 

CO 

1, 

t — 1 

CO 

81, 

t — 1 

CO 

81 

$BCINP  ISPARTOF  = 

'BLK- 

3', 

IBTYP  = 

5, 

5, 

5, 

40, 

40, 

5 

IBDIR  = 

2, 

3, 

-3, 

-1, 

1, 

-2 

JBCS  = 

1, 

1, 

1, 

61, 

1, 

1 

JBCE  = 

61, 

61, 

61, 

61, 

1, 

61 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

81 

KBCE  = 

1, 

81, 

81, 

t — 1 

CO 

CO 

81 

LBCS  = 

1, 

1, 

81, 

1, 

1, 

1 

LBCE  = 
$END 

t — 1 

CO 

1, 

81, 

t — 1 

CO 

CO 

81 

$BCINP  ISPARTOF  = 

'BLK- 

■4', 

IBTYP  = 

5, 

5, 

5, 

40, 

5, 

5 

IBDIR  = 

2, 

-1, 

3, 

1, 

-3, 

-2 

JBCS  = 

1, 

61, 

1, 

1, 

1, 

1 

JBCE  = 

t — 1 

61, 

t — 1 

1, 

t — 1 

61 

90 


KBCS  = 

1, 

1, 

1, 

1, 

1, 

81 

KBCE  = 

1, 

t — 1 

CO 

t — 1 

CO 

CO 

81, 

81 

LBCS  = 

1, 

1, 

1, 

1, 

81, 

1 

LBCE  = 

CO 

CO 

1, 

CO 

81, 

81 

$END 

$BCINP  ISPARTOF  =  'BLK-5', 

IBTYP  =  5, 

IBDIR  =  2, 

JBCS  =  1, 

JBCE  =  41, 

KBCS  =  1, 

KBCE  =  1, 

LBCS  =  1, 

LBCE  =  31, 

$END 

$BCINP  ISPARTOF  =  'BLK-6', 

IBTYP  =  5, 

IBDIR  =  2, 

JBCS  =  1, 

JBCE  =  41, 

KBCS  =  1, 

KBCE  =  1, 

LBCS  =  1, 

LBCE  =  31, 

$END 


$BCINP  ISPARTOF  =  'BLK-7', 


IBTYP  = 

5, 

40, 

40, 

40 

IBDIR  = 

2, 

1, 

3, 

-3 

JBCS  = 

1, 

1, 

1, 

1 

JBCE  = 

45, 

1, 

45, 

45 

KBCS  = 

1, 

1, 

1, 

1 

KBCE  = 

1, 

t — 1 

61, 

61 

LBCS  = 

1, 

1, 

1, 

31 

LBCE  = 
$END 

\ — 1 

00 

31, 

1, 

31 

$BCINP  ISPARTOF  = 

'BLK 

-8', 

IBTYP  = 

5, 

40, 

40, 

40 

IBDIR  = 

2, 

-1, 

3, 

-3 

JBCS  = 

1, 

45, 

1, 

1 

JBCE  = 

45, 

45, 

45, 

45 

KBCS  = 

1, 

1, 

1, 

1 

KBCE  = 

1, 

61, 

t — 1 

61 

LBCS  = 

1, 

1, 

1, 

31 

LBCE  = 
$END 

\ — 1 

00 

31, 

1, 

31 

$BCINP  ISPARTOF  = 

'BLK 

-9', 

IBTYP  = 

5, 

40, 

40, 

40 

IBDIR  = 

2, 

-1, 

-3, 

3 

JBCS  = 

1, 

45, 

1, 

1 

JBCE  = 

45, 

45, 

45, 

45 

KBCS  = 

1, 

1, 

1, 

1 

KBCE  = 

1, 

61, 

61, 

61 

LBCS  = 

1, 

1, 

31, 

1 

LBCE  = 
$END 

31, 

31, 

31, 

1 

$BCINP  ISPARTOF  = 

'BLK- 

10  '  , 

IBTYP  = 

5, 

40, 

40, 

40 

IBDIR  = 

2, 

1, 

-3, 

3 

91 


JBCS 

= 

1, 

1,  1, 

1, 

JBCE 

45, 

1,  45, 

45, 

KBCS 

1, 

1,  1, 

1, 

KBCE 

1, 

61,  61, 

61, 

LBCS 

1, 

1,  31, 

1, 

LBCE 

$END 

31, 

31,  31, 

1, 

$BCINP  ISPARTOF  = 

'BLK-11 ' , 

IBTYP 

= 

5, 

40,  40, 

40, 

40, 

IBDIR 

2, 

-3,  -3, 

1, 

-1, 

JBCS 

1, 

1,  25, 

1, 

41, 

JBCE 

41, 

25,  41, 

1, 

41, 

KBCS 

1, 

1,  1, 

1, 

1, 

KBCE 

1, 

61,  61, 

61, 

61, 

LBCS 

1, 

45,  45, 

1, 

1, 

LBCE 

$END 

45, 

45,  45, 

45, 

45, 

$BCINP  ISPARTOF  = 

'BLK-12 ' , 

IBTYP 

= 

5, 

40,  40, 

40, 

40, 

IBDIR 

2, 

-3,  -3, 

-1, 

1, 

JBCS 

1, 

17,  1, 

41, 

1, 

JBCE 

41, 

41,  17, 

41, 

1, 

KBCS 

1, 

1,  1, 

1, 

1, 

KBCE 

1, 

61,  61, 

61, 

61, 

LBCS 

1, 

45,  45, 

1, 

1, 

LBCE 

$END 

45, 

45,  45, 

45, 

45, 

$BCINP  ISPARTOF  = 

'BLK-13 ' , 

IBTYP 

= 

5, 

5,  40, 

40, 

40, 

40 

IBDIR 

2, 

2,  -1, 

-3, 

1, 

3 

JBCS 

1, 

1,  45, 

1, 

1, 

1 

JBCE 

45, 

45,  45, 

45, 

1, 

45 

KBCS 

1, 

1,  1, 

1, 

1, 

1 

KBCE 

1, 

1,  61, 

61, 

61, 

61 

LBCS 

51, 

1,  1, 

65, 

1, 

1 

LBCE 

$END 

65, 

51,  65, 

65, 

65, 

1 

$BCINP  ISPARTOF  = 

'BLK-14 ' , 

IBTYP 

= 

5, 

5,  40, 

40, 

40, 

IBDIR 

2, 

2,  -1, 

1, 

3, 

JBCS 

1, 

1,  41, 

1, 

1, 

JBCE 

41, 

41,  41, 

1, 

41, 

KBCS 

1, 

1,  1, 

1, 

1, 

KBCE 

1, 

1,  61, 

61, 

61, 

LBCS 

51, 

1,  1, 

1, 

1, 

LBCE 

$END 

65, 

51,  65, 

65, 

1, 

$BCINP  ISPARTOF  = 

'BLK-15 ' , 

IBTYP 

= 

5, 

5,  40, 

40, 

40, 

40 

IBDIR 

2, 

2,  -1, 

-3, 

1, 

3 

JBCS 

1, 

1,  45, 

1, 

1, 

1 

JBCE 

45, 

45,  45, 

45, 

1, 

45 

KBCS 

1, 

1,  1, 

1, 

1, 

1 

KBCE 

1, 

1,  61, 

61, 

61, 

61 

LBCS 

51, 

1,  1, 

65, 

1, 

1 

LBCE 

$END 

= 

65, 

51,  65, 

65, 

65, 

1 

$BCINP  ISPARTOF  =  'BLK-16', 
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IBTYP  = 

40, 

5, 

5, 

40, 

40, 

40 

IBDIR  = 

3, 

2, 

2, 

1, 

-1, 

-3 

JBCS  = 

1, 

1, 

1, 

1, 

45, 

1 

JBCE  = 

45, 

45, 

45, 

1, 

45, 

45 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

1 

KBCE  = 

t — 1 

1, 

1, 

t — 1 

t — 1 

61 

LBCS  = 

1, 

51, 

1, 

1, 

1, 

65 

LBCE  = 
$END 

1, 

65, 

t — 1 

LO 

65, 

65, 

65 

$BCINP  ISPARTOF  = 

'BLK- 

■17  '  , 

IBTYP  = 

40, 

5, 

5, 

40, 

40, 

IBDIR  = 

3, 

2, 

2, 

1, 

-1, 

JBCS  = 

1, 

1, 

1, 

1, 

41, 

JBCE  = 

41, 

41, 

41, 

1, 

41, 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

KBCE  = 

t — 1 

1, 

1, 

61, 

t — 1 

LBCS  = 

1, 

1, 

51, 

1, 

1, 

LBCE  = 
$END 

1, 

t — 1 

LO 

65, 

65, 

65, 

$BCINP  ISPARTOF  = 

'BLK- 

CO 

IBTYP  = 

40, 

5, 

5, 

40, 

40, 

40 

IBDIR  = 

3, 

2, 

2, 

1, 

-1, 

-3 

JBCS  = 

1, 

1, 

1, 

1, 

45, 

1 

JBCE  = 

45, 

45, 

45, 

1, 

45, 

45 

KBCS  = 

1, 

1, 

1, 

1, 

1, 

1 

KBCE  = 

t — 1 

1, 

1, 

61, 

61, 

61 

LBCS  = 

1, 

1, 

51, 

1, 

1, 

65 

LBCE  = 
$END 

1, 

t — 1 

LO 

65, 

65, 

65, 

65 

$BCINP  ISPARTOF  =  'BLK-19', 

IBTYP  =  5, 

IBDIR  =  2, 

JBCS  =  1, 

JBCE  =  -1, 

KBCS  =  1, 

KBCE  =  1, 

LBCS  =  1, 

LBCE  =  -1, 

$END 

B.3  Boundary  condition  file 


52 

3 

1  61 
1  81 
1  1 
1 
1 
-3 

1  61 
1  81 
81  81 
92 
-1 

61  61 
1  81 
1  81 
1 
2 


number  of  boundary  surfaces  ===>  BLK#-1 
!***  periodic  MB  #1=  type  #52 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  patched  MB  #2=  type  #92 

!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 
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1 

1 

1 

1 

52 

-3 

1 

1 

81 

1 

1 

3 

1 

1 

1 

22 

2 

1 

1 

1 

0 

14 

1 

1 

1 

1 

0 

13 

-2 

1 

81 

1 

0 

6 

22 

2 

1 

1 

1 

0 

41 

3 

1 

1 

1 

0 

92 

-1 


1 

3 

1 

1 

1 

1 

41 

-3 

1 

1 

81 

0 

92 


1 

81 

81 


61 

81 

81 


61 

81 

1 


61 

1 

81 


1 

81 

81 


61 

81 

81 


61 

1 

81 


61 

81 

1 


1 

81 

81 


61 

81 

81 


coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  periodic  MB  #1=  type  #52 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Prescribed  =  type  #14 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Farfield  (dp/dn=0)  =  type  #13 

coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
number  of  boundary  surfaces  ===>  BLK#-2 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Periodic  =  type  #41 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #3=  type  #92 
coordinate  direction  normal  to  surface 


61 

61 

1 

start , 

end 

in 

i-direction 

1 

81 

1 

start , 

end 

in 

j-direction 

1 

81 

1 

start , 

end 

in 

k-direction 

flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  Periodic  =  type  #41 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #1=  type  #92 
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1 

1  1 
1  81 
1  81 
1 
1 

-1 

61  61 
1  81 
1  81 
13 
-2 

1  61 
81  81 
1  81 
0 
6 

22 

2 

1  61 
1  1 

1  81 
0 
52 
3 

1  61 
1  81 
1  1 

1 
3 
-3 

1  61 
1  81 
81  81 
52 
-3 

1  61 
1  81 
81  81 
1 
3 

3 

1  61 
1  81 
1  1 
92 
-1 

61  61 
1  81 
1  81 
1 

4 
1 

1  1 

1  81 
1  81 
92 
1 

1  1 

1  81 
1  81 
1 
2 

-1 


coordinate  direction  normal  to  surface 

start,  end  in  i-direction 

start,  end  in  j-direction 

start,  end  in  k-direction 

flag  for  order  of  first  derivative  be 

donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  Farfield  (dp/dn=0)  =  type  #13 

coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
number  of  boundary  surfaces  ===>  BLK#-3 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  periodic  MB  #3=  type  #52 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  periodic  MB  #3=  type  #52 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #4=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #2=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
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61 

61 

1 

81 

1 

81 

13 

-2 

1 

61 

81 

81 

1 

81 

0 

6 

22 

2 

1 

61 

1 

1 

1 

81 

0 

11 

-1 

61 

61 

1 

81 

1 

81 

0 

41 

3 

1 

61 

1 

81 

1 

1 

0 

92 

1 

1 

1 

1 

81 

1 

81 

1 

3 

-1 

61 

61 

1 

81 

1 

81 

41 

-3 

1 

61 

1 

81 

81 

81 

0 

13 

-2 

1 

61 

81 

81 

1 

81 

0 

1 

22 

2 

1 

41 

1 

1 

1 

31 

0 

1 

22 

2 

1 

41 

1 

1 

1 

31 

!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  Farfield  (dp/dn=0)  =  type  #13 

!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  number  of  boundary  surfaces  ===>  BLK#-4 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  Exit  =  type  #11 

!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  Periodic  =  type  #41 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  patched  MB  #3=  type  #92 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  Periodic  =  type  #41 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  Farfield  (dp/dn=0)  =  type  #13 

!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  number  of  boundary  surfaces  ===>  BLK#-5 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  number  of  boundary  surfaces  ===>  BLK#-6 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 
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flag  for  order  of  first  derivative  be 
number  of  boundary  surfaces  ===>  BLK#-7 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #8=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #16=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #11=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-8 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #7=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #15=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
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number  of  boundary  surfaces  ===>  BLK#-9 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #10=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #11=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #18=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-10 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
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***  patched  MB  #9=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #12=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #13=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-11 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #12=  type  #92 
coordinate  direction  normal  to  surface 
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flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #12=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #7=  type  #92 
coordinate  direction  normal  to  surface 
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start,  end  in  i-direction 

start,  end  in  j-direction 

start,  end  in  k-direction 

flag  for  order  of  first  derivative  be 

donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #9=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-12 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #11=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #11=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #8=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #10=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
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start,  end  in  k-direction 

flag  for  order  of  first  derivative  be 

donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-13 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #14=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #10=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #18=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
number  of  boundary  surfaces  ===>  BLK#-14 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
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start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #15=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #13=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
number  of  boundary  surfaces  ===>  BLK#-15 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #16=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #8=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
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!  start,  end  in  k-direction 
!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  patched  MB  #14=  type  #92 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  number  of  boundary  surfaces  ===>  BLK#-16 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  patched  MB  #15=  type  #92 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  patched  MB  #17=  type  #92 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 
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1 

65 

!  donor  mesh 

start,  end  in  k-direction 

92 

! ***  patched 

MB  #7=  type  #92 

-3 

!  coordinate 

direction  normal  to  surface 

1 

45 

!  start,  end 

in  i-direction 

1 

61 

!  start,  end 

in  j-direction 

65 

65 

!  start,  end 

in  k-direction 

1 

!  flag  for  order  of  first  derivative  be 

7 

!  donor  mesh 

3 

!  coordinate 

direction  normal  to  surface 

1 

45 

!  donor  mesh 

start,  end  in  i-direction 

1 

61 

!  donor  mesh 

start,  end  in  j-direction 

1 

1 

!  donor  mesh 

start,  end  in  k-direction 

5 

22 

3 

1 

1 

1 

0 

22 

2 

1 

1 

1 

0 

22 

2 

1 

1 

51 

0 

92 

1 

1 

1 

1 

1 

16 

-1 

45 

1 

1 

92 

-1 

41 

1 

1 

1 

18 

1 

1 

1 

1 

6 

22 

3 

1 

1 

1 

0 

22 

2 

1 

1 


41 

61 

1 


41 

1 

51 


41 

1 

65 


1 

61 

65 


45 

61 

65 


41 

61 

65 


1 

61 

65 


45 

61 

1 


45 

1 


number  of  boundary  surfaces  ===>  BLK#-17 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  patched  MB  #16=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
***  patched  MB  #18=  type  #92 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
donor  mesh 

coordinate  direction  normal  to  surface 
donor  mesh  start,  end  in  i-direction 
donor  mesh  start,  end  in  j-direction 
donor  mesh  start,  end  in  k-direction 
number  of  boundary  surfaces  ===>  BLK#-18 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
start,  end  in  k-direction 
flag  for  order  of  first  derivative  be 
***  Rel.  Frame  No-slip  =  type  #22 
coordinate  direction  normal  to  surface 
start,  end  in  i-direction 
start,  end  in  j-direction 
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51 


1 

0 

22 

2 

1 

1 

51 

0 

92 

1 

1 

1 

1 

1 

17 

-1 

41 

1 

1 

92 

-1 

45 

1 

1 

1 

13 

1 

1 

1 

1 

92 

-3 

1 

1 

65 

1 

9 

3 

1 

1 

1 


45 

1 

65 


1 

61 

65 


41 

61 

65 


45 

61 

65 


1 

61 

65 


45 

61 

65 


45 

61 

1 


!  start,  end  in  k-direction 
!  flag  for  order  of  first  derivative  be 
!***  Rel.  Frame  No-slip  =  type  #22 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!***  patched  MB  #17=  type  #92 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  patched  MB  #13=  type  #92 
!  coordinate  direction  normal  to  surface 
!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 

!***  patched  MB  #9=  type  #92 
!  coordinate  direction  normal  to  surface 

!  start,  end  in  i-direction 

!  start,  end  in  j-direction 

!  start,  end  in  k-direction 

!  flag  for  order  of  first  derivative  be 
!  donor  mesh 

!  coordinate  direction  normal  to  surface 
!  donor  mesh  start,  end  in  i-direction 

!  donor  mesh  start,  end  in  j-direction 

!  donor  mesh  start,  end  in  k-direction 
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